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I.  INTRODUCTION 

This  report  is  designed  to  give  the  reader  an  in-depth  view  of  the  development  of  an 
adaptive.  Kalman,  target  tracking  filter.  Emphasis  is  placed  on  the  synthesis  and  analysis  so 
as  to  convey  the  rationale  for  various  design  decisions.  Since  this  filter  was  specifically  de- 
veloped for  the  improved  (digital)  MARK  68  Gunfire  Control  System  (GFCS).  the  parametric 
studies  necessarily  teflect  values  characteristic  to  such  systems  and  their  particular  scenarios.* 
To  the  extent  possible,  however,  the  analysts  has  been  generalized  and  deals  in  a broad 
manner  with  the  problem  of  tracking  and  predicting  the  stale  of  a general  target  for  fire 
control  or  other  purposes.  The  MARK  68  GFCS  is  in  fact  multifunctional  as  it  must  be 
able  to  track  and  engage  targets  of  ail  types.  For  the  most  part,  the  various  techniques  used 
in  this  filter/predictor  are  not  new  and  can  be  found  scattered  througnout  the  literature,  as 
can  be  seen  from  the  rather  lengthy  list  of  references.  It  is  hoped,  however,  that  this  report 
will  be  a useful  reference  for  others  who  will  work  in  this  area  in  that  it  combines  these 
various  state-of-the-art  concepts  and  brings  them  to  bear  on  this  particular  application.  Tins 
document  summarizes  work  performed  to  date  and  also  indicates  the  directions  that  future 
investigations  might  take. 


Let  us  first  consider  the  bask:  functions  performed  bv  a gunfire  control  system.  As 
shown  in  Figure  1.1.  the  heart  of  the  system  is  the  so-called  predictor-ballistics  loop.  Es- 
sentially. the  ballistics  section  generates  the  possible  four-dimensional  trajectories  (or 


Figure  1.1.  Simple  Functional  Diagram  of  a Gunfire  Control  System 
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perhaps  only  terminal  coordinates)  of  ownship’s  projectiles.  The  prediction  section  must 
supply  sintilar  information  about  the  target  trajectory  . The  loop  is  then  iteratively  closed 
by  matching  the  spatial  and  temporal  coordinates  of  the  projectile  to  the  target  to  effect  an 
intercept.  In  this  report,  we  will  be  dealing  only  with  the  ’“front  end"  of  the  system,  i.e.. 
from  the  target  to  the  predictor.  Quite  obviously  the  predictor  is  the  ultimate  product  re- 
quired for  the  front  end  of  the  loop,  and  this  factor  will  be  repeatedly  stressed  in  this  report. 
If  we  had  some  other  means  of  specifying  the  target  trajectory,  there  would  be  no  need  for 
the  other  elements  of  the  front  end.  In  fact,  of  course  the  target  trajectory  is  not  directly 
available  to  the  Fire  control  system.  Instead,  the  sensor  tracks  the  target’s  current  position, 
superimposing  measurement  error  in  so  doing.  The  purpose  of  the  filter  then  is  to  process 
the  noisy  position  measurements  in  such  a manner  as  to  estimate  the  parameters  required 
for  the  prediction  model.  Such  parameters  might  include  smooth  current  Cargct  position, 
velocity,  and  acceleration. 

it  is  assumed  for  this  report  that  measurements  are  made  of  target  position  only  since 
rates  (e.g.  Doppler,  gyros,  etc.)  are  not  currently  available  from  GFCS  sensors.  I?  will  Ire 
assumed  that  these  measurements  have  been  stabilized  (i.e..  ownship  angular  motion  removed) 
and  transformed  to  Cartesian  coordinates.  The  implications  of  this  transformation  will  be 
considered  at  length  latet  in  the  report.  Tire  use  of  stabilized  coordinates,  however,  is 
probably  universal  since  one  does  not  want  the  target  motion  filter  to  have  to  deal  with 
ownship  angular  motion.  The  origin  of  the  coordinate  system  shall  be  a reference  point  on 
ownship  and  common  to  all  sensors.  Tracking,  filtering,  and  prediction  wili  therefore  be 
performed  in  ownship  coordinates.  This  "ownship  coordinate  system'’  is  defined  simply  as 

Xt  Jil  * Xr(t)  - Xodt)  (1.1) 


where  Xt <*)  and  X»,<U  are  the  motion  of  the  target  and  of  ownship  relative  to  some  earth- 
fixed  or  inertial  reference  frame.'’  Whenever  it  is  necessary  to  transform  to  an  inertial  frame, 
say  for  purposes  of  ballistic  computations,  or  to  account  (if  necessary)  for  ownship  maneu- 
vers, Equation  (1 . 1 1 and  its  derivatives  would  bo  used.  The  reason  for  this  choice  of 
Cartesian  ownship  coordinates  for  target  filtering  and  prediction  b simply  that  the  motion 
of  most  targets  in  a GFCS  scenario  is  belter  behaved  or  more  closely  linear-in  this  system 
than  ir,  others.  One  must  consider  the  fact  that,  if  a hostile  target  b so  close  as  to  be  in  the 
range  of  fire  of  a gun  system,  most  likely  the  engagement  is  mutual  in  that  both  ownship 
and  target  are  maneuvering  retatir<  to  each  other  in  order  to  effect  their  own  fire  control 
strategy. 

Let  us  now  consider  lire  nature  of  this  problem  from  the  point  of  view  of  error  sources 
and  constraints  on  the  solution  method.  The  prediction  error,  that  is  the  error  in  determin- 
ing the  projectile-target  intercept  point,  basically  results  from  two  (not  necessarily  independent) 
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sources:  predictor  modeling  error  and  filtering  or  estimation  error.  Prediction  modeling 
error  occurs  because,  not  knowing  the  target  strategy,  an  incorrect  functional  form  of  the 
predictor  is  utilized.  While  we  strive  to  construct  a reasonable  approximation  for  the  actual 
prediction  model  of  the  general  unknown  target,  we  can  never  do  this  exactly.  Therefore, 
even  if  we  effected  perfect  filtering  (i.e.,  perfect  estimation  of  target  parameters  such  as  cur- 
rent position  and  rates),  we  would  nevertheless  suffer  prediction  errors  due  to  unaccountable 
target  maneuvers  during  the  prediction  time.  Of  course,  even  if  target  modeling  were  perfect, 
there  would  remain  estimation  errors  due  to  sensor  measurement  noise  which  would  be 
extrapolated  from  current  time  by  the  predictor.  The  dominant  variable,  an  effect  which 
significantly  determines  the  relative  and  absolute  magnitude  of  these  prediction  errors,  is  the 
prediction  time.  We  shall  see  that  prediction  errors  are  magnified  exponentially  with  in- 
creased prediction  time. 

Target  modeling  error  contributes  to  the  filter  error  as  well  as  the  prediction  error, 
particularly  for  maneuvering  targets.  In  general,  a target  will  probably  be  maneuvering  in 
some  unknown  fashion-at  least  as  far  as  the  GFCS  is  concerned  -when  the  target  is  in  range 
of  ownship  guns.  There  are  several  reasons  for  this.  First,  the  target  probably  has  a mission 
to  accomplish  or  it  would  not  be  there.  The  target  has  a particular  strategy  to  achieve  that 
goal:  and  very  often  that  goal  and  strategy  are  known  to  members  of  ownship  crew.  Due  to 
the  complexity  of  this  information,  however,  such  knowledge  is  not  used  (except  very 
indirectly)  by  current  gun-fire  control  systems.  Examples  of  such  targets  are  anti-ship 
missiles  that  terminally  home  on  their  targets  or  other  weapons  platforms  (such  as  fighter 
bombers),  which  must  maneuver  in  specific  ways  in  order  to  solve  their  fire  control  problem. 
(The  goal-oriented  target  and  its  possible  prediction  by  a fire  control  system  are  the  subject 
of  another  report  to  be  published.)  Another  important  reason  for  target  maneuvers  is  to 
evade  ownship  weapons  Certainly,  manned  targets  in  range  of  our  guns  would  bt  well 
aware  of  the  principles  of  GFC'S  operation  and  behave  accordingly.  It  is  also  probable  that 
missiles  will  be  made  capable  of  evasive  action  either  by  preprogramming  or  by  being  under 
remote  control.  There  are  other  reasons  for  target  maneuver  such  as  atmospheric  turbulence, 
target  sensor-induced  errors,  target  control  system  behavior  etc.  In  any  case,  it  is  not  pos- 
sible to  fully  and  deterministically  model  the  anticipated  motion  of  a general  target,  and  we 
must  live  with  a certain  amount  of  tracking  and  prediction  error  due  to  this  limitation.  We 
shall  find  that  there  are  fundamental  limitations  on  the  ability  to  track  maneuvering  targets. 

A maneuver  will  subsequently  be  defined  for  purposes  of  this  report  as  any  target  motion 
that  tends  to  cause  a certain  filter  implementation  to  diverge  (i.e..  estimation  errors  become 
large). 

Several  other  factors  affect  the  filter  estimation  error  and  these  must  each  he  considered. 
Tracking  sensor  accuracy  is  an  important  one.  A very  important  matter  is  filter  settling 
time.  In  a situation  with  more  than  one  target  to  be  engaged  or  when  target  acquisition  and 
tracking  commence  late,  it  is  mandatory  that  the  filter  settle  quickly,  within  a few  seconds. 
Computational  constraints  on  such  factors  as  data  rate  and  model  complexity  must  be 


considered.  Computer  word  length  and  real-time  implementation  constraints  pose  severe 
limitations  on  what  one  might  hope  to  accomplish  when  dealing  with  on-line  filtering  and 
prediction.  In  short,  the  problem  of  target  estimation  for  gunfire  control  application  must 
be  considered  in  toto,  and  in  so  doing,  the  filter  designer  finds  himself  facing  many  limita- 
tions, both  fundamental  and  practical,  that  must  be  recognized  if  a reasonable  design  is  to 
be  realized.  In  thf*  final  analysis,  judgments  must  be  made  that,  on  occasions,  are  difficult 
to  justify  with  theoretical  rigor,  and  filter  design  becomes  somewhat  an  art  as  well  as  a 
science. 

Before  commencing  the  actual  report,  a brief  discussion  of  the  nature  and  structure  of 
this  study  is  in  order.  The  investigation,  to  the  extent  possible,  is  primarily  parametric  in 
that  the  effects  of  individual  parameters  affecting  performance  were  isolated,  wherever  pos- 
sible. in  order  to  assess  the  sensitivity  of  performance  to  variation  of  the  parameters.  It 
should  be  emphasized  that  many  of  the  results  presented  in  this  report  were  based  on  studies 
designed  to  reveal  the  relative  merits  of  one  technique  or  parameter  with  respect  to  another 
and,  as  such,  ma--  lot  always  be  indicative  of  actual  implemented  performance  when  in- 
corporated wit.i  oinri  techniques  and  parameters  in  an  actual  tracking  situation.  The 
criterion  for  performance,  or  “norm.”  was  chosen,  due  to  the  nature  of  the  gun  fire  control 
problem,  to  be  prediction  performance.  For  a “nominal”  Navy  gun,  a time  of  flight  (or 
prediction  time)  of  1 0 seconds  was  chosen  as  typical  for  engagement  of  air  targets  and  is 
used  consistently  throughout  this  report.  The  actual  statistical  form  of  the  prediction  error, 
used  to  evaluate  the  relative  merits  of  parameters  or  filter  formulations,  depends  upon  the 
situation.  Whenever  practical,  the  error  covariance  approach  was  utilized  since  it  essentially 
represents  the  result  of  a perfect  Monte  Carlo  study  (repeated  simulation  trials  plus 
averaging).  The  error  covariance  approach  can  be  used  to  evaluate  optimal  or  suboptimal  (if 
the  “real  world"  or  truth  model  is  known)  performance.  There  are  many  cases,  however, 
when  the  truth  model  is  not  readily  available,  and  simulation  must  be  employed  with  root- 
mean-square  prediction  errors  used  as  a norm.  This  is  particularly  true  in  evaluating  a non- 
linear adaptive  filter  against  various  nonlinear  trajectories.  Ultimately,  of  course,  the  filter 
must  be  evaluated  in  relatively  uncontrolled  circumstances  using  actual  tracking  data  of  real 
target  trajectories. 

Finally,  let  us  consider  the  spatial  problem.  Measurements  of  target  position  arc  made 
in  a spherical  polar  coordinate  frame,  but  targets  are  generally  not  linearly  described  in  such 
a frame.  Therefore,  a nonlinearity  is  going  to  appear  in  the  full  three-dimension  problem. 
This  consideration  was  postponed  until  the  last  chapter,  and  all  sections  until  then  discuss 
only  one-dimensional  estimation  of  target  motion  as  a function  of  time.  The  reader  may 
consider  this  single  “channel"  to  be  one  direction  of  some  orthogonal  system  to  be  defined 
later. 

Since  the  author  has  noticed  an  occasional  tendency  among  persons  unfamilar  with 
Kalman  filtering  to  regard  the  technique  as  exceedingly  complex,  mysterious  and/or 
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omnipotent,  he  feels  compelled  to  advise  the  reader  that  the  antithesis  is  probably  much 
closer  to  the  truth.  In  the  next  section,  the  Kalman  filter  will  be  introduced  and,  throughout 
this  report,  the  reader  will  note  that  performance  remains  limited,  as  it  must  be,  by  the 
quality  and  quantity  of  information  available. 


II.  GENERAL  DESCRIPTION  OF  THE  DISCRETE  KALMAN  FILTER 

In  this  section,  the  concepts  and  equations  of  Kalman  filters  will  be  very  briefly  re- 
viewed in  order  to  familiarize  the  reader  with  the  notation  and  the  types  of  information  re- 
quired to  construct  a filter  of  this  type.  The  theoretical  elegance  and  practical  utility  of  this 
estimator  have  made  it  a popular  subject  for  both  mathematicians  and  engineers  in  recent 
years.  For  the  reader  who  wishes  to  pursue  the  subject  in  greater  depth,  the  author  recom- 
mends such  literature  as  Gelb  ( 1 974),  Sorenson  ( 1 966)  or  one  of  many  other  textbooks  that 
deal  with  estimation  theory  such  as  Morrison  (1969)  or  Sage  and  Melsa  (1971).*  Wc  will 
deal  here  only  with  the  discrete  formulations  since  they  arc  the  more  natural  for  digital 
computer  implementation. 

Before  considering  the  equations,  we  might  first  qualitatively  define  a Kalman  filter. 

The  Kalman  filter  is  a linear  minimum-variance  estimator.  In  fact,  it  is  a recursive  formula- 
tion of  the  minimum  variance  rules  for  combining  certain  a priori  information  with  a sequence 
of  observations  that  contain  noise.  We  might  also  point  out  now  that  the  required  a priori 
information  must  usually  be  chosen  with  considerable  care  if  one  is  to  achieve  the  desired 
results.  A Kalman  filter  is  also  a maximum-likelihood  estimator  (and  in  that  sense  an  optimal 
linear  or  nonlinear  filter)  but  only  when  the  errors  are  uncorrelated  and  Gaussian.  Of 
course,  as  mentioned  in  the  Introduction,  we  really  do  not  expect  to  operate  our  filter  under 
optimal  conditions  due  primarily  to  target  modeling  errors  and  also  to  correlated  and  per- 
haps slightly  non-Gaussian  measurement  errors.  Indeed,  most  applications  of  the  Kalman 
filter  are  suboptimal.  Kalman  filters  often  are  actually  designed  to  be  suboptimal.  This 
frequently  occurs  when  the  dimension  of  the  true  state  vector  is  so  large  as  to  render  the 
computations  impossible  to  effect  on  a selected  digital  computer.  Of  course,  we  strive  to 
approach  optimality  within  the  imposed  computational  constraints  and  do  just  that  when 
circumstances  permit.  The  fact  that  wc  can  often  approach  optimal  filter  performance  is 
therefore  one  important  reason  for  selecting  a Kalman  filter. 

Another  equally  important  reason  is  the  inherent  flexibility  and  control  one  has  over 
the  filtering  process  in  order  to  maintain  desired  or  expected  filter  performance.  A signifi- 
cant computational  advantage  over  certain  other  filtering  techniques  (such  as  a Baycscan 
estimator)  is  also  enjoyed  by  Kalman  filters  due  to  a reduced  matrix  inversion  burden.  Singer 
and  Behnke  (1971 ),  in  one  of  the  best  papers  on  the  subject,  compare  five  types  of  filters 


‘See  References. 
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in  terms  of  tracking  accuracy  and  computational  requirements  when  tracking  maneuvering  I 

targets  in  a tactical  application.  They  considered  first  and  second  order  polynomial  Kalman  1 

filters  (to  be  discussed  in  the  next  section);  an  alpha-beta  filter  (a  recursive  least  squares  * I 

filter  with  a first  order  target  model);  a Wiener  filter  (steady-state  Kalman  filter);  and  a 

simple  two-point  extrapolator.  The  second-order  Kalman  filter  was  found  to  be  the  most 

accurate.  In  a later  section,  we  shall  also  demonstrate  that  it  is  possible  to  significantly 

reduce  by  more  than  an  order  of  magnitude -the  computational  burden  of  a Kalman  filter. 

It  is  generally  accepted  that,  under  optimal  conditions,  a Kalman  filter  will  out-perform 
other  data-processing  algorithms  such  as  those  based  upon  least  squares.  This  factor,  although 
important,  is  probably  not  in  itself  justification  for  implementing  a Kalman  filter  since, 
under  optimal  conditions,  almost  any  algorithm  will  perform  acceptably.  The  real  advan- 
tage of  a Kalman  filter  is  under  suboptimal  conditions.  The  Kalman  filter  “thinks"  it  knows 
how  well  it  is  performing  and  can  be  made  to  determine  when  it  is  not  performing  properly. 

Then,  as  just  mentioned,  the  Kalman  filter  can  utilize  its  flexibility  to  adjust  the  covariance 
and  regain  calculated  performance.  These  concepts  will  be  discussed  in  much  greater  depth 
in  the  section  on  adaption. 

The  target  is  assumed  to  be  modeled  by  a linear,  discrete,  system  model,  as  shown  in 
Figure  2.1(a).  The  equation  for  the  system  is: 

xfk  + I)  = p(k  + l.  k)x(k)  + w(k)  (2.1) 

The  system  at  time  t(k)  is  characterized  by  a state  vector  xfk)  of  dimension  n.  The  state 
vector  is  assumed  to  p-opagate  linearly  according  to  an  (n  X n)  transition  matrix  0(k  + I.  k). 

For  purposes  of  this  presentation,  no  control  input  or  known  forcing  function  will  be 
considered. 

Any  unmodeled  effects  will  be  assumed  to  be  random  or  "nondetcrministic"  and  are 
lumped  into  a term  (the  last  tertn-w)  called  “process  noise."  The  random  process  noise 
sequence  has  known  statistics  given  by 

H |w(k)!  = 0 (2.2) 

F(w(j)w'(k)|  = Q(k)6ji,  (2.3)  < 

i 

The  actual  magnitude  of  Q is  required  by  the  Kalman  filter  for  proper  (optimal)  operation. 

We  shall  see  that  often  this  information  is  not  available  a priori  and  0 must  be  adjusted  dur- 
ing operation. 

In  Figure  2.1(b).  we  assume  we  have  measurements  ..(k)  of  dimension  m,  related  to  the 
state  vector  by  an  (m  X n)  “observation"  matrix  H(k)  (also  called  the  “matrix  of  partial 
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derivatives")  corrupted  with  a zero-mean.  Gaussian,  white  noise  sequence  v(k).  That  is. 
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z(k)  = H(k)  x(k)  + v(k)  (2.4) 

where 

E[v(k)]  = 0 (2.5) 

Elv(j)vT(k)  = R(k)6jk  (2.6) 

(Notice  that  the  use  of  the  H matrix  in  Equation  (2.4)  readily  allow.-  the  mixing  of  measure- 
ments of  different  types.] 

It  is  further  assumed  that 


E(v(j)wT(k)]  = 0 (2.7) 

The  basic  recursive,  linear  estimation  problem  is  to  determine  an  estimate  x(kjk)  of 
x(k)  that  is  a linear  combination  of  the  previous  estimate  and  the  current  measurement  as 
in  Figure  2.1(cv 


We  find  we  can  write  this  as 


x(kik)  = xfklk-  1)  + K(k) i^kik - !)  (2.8) 

where 

x(kik-  1)  = 0(k.  k - I)  xCk-  lik  - l)  (2.9) 

is  the  extrapolated  estimate  at  t(k).  i.e..  the  estimate  at  t(k)  based  on  measurement  data 
through  t(k  - 1 ).  The  a ptl'tri  residual  vector  of  dimension  m is  defined 


} £(kik-  I)  = z(k)  - H(k)x(k|k-  1)  (2.10 

» 

I and  k(k>  is  the  (n  X in)  gain  matrix.  The  difference  between  the  estimate  and  the  true  state 

| is  the  error  vector 
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e(jik)  = x(jlk)  - x(j)  (2.11) 

with  associated  (n  X nl  error  covariance  matrix  defined  as 

P(jik)  = E(c(jik)  eT(jik)|  (2.«2) 

The  basis  of  ’he  Kalman  filter  is  the  selection  of  the  gain  matrix  in  Equation  (2.8)  that 
minimizes  the  trace  of  the  error  covariance  matrix.  The  solution  to  this  problem,  found  by 
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Kalman  (I960.  1961 ).  is  well  known  and  its  derivation  can  be  found  in  the  previously  men- 
tioned references.  The  Kalman  filter  is  found  to  be  a set  of  recursive  matrix  equations  thaj 
are  usually  written  and  implemented  in  two  distinct  steps.  The  "extrapolation”  process 
simply  predicts  the  state  vector,  using  Equation  (2.9),  and  the  error  covariance  via 


Pfkik-  1)  = #k.k-  1)P(k-  Ilk-  I) 0T(k, k - I)  + Q(k  - l; 


(2. 13 1 


in  this  step,  the  error  covariance  usually  increases  since  we  are  advancing  in  time  without 
benefit  of  new  measurements.  The  second  step,  the  “update."  applies  to  the  process  of 
recstimating  the  system  to  acknowledge  that  new  information.  z(k).  is  available.  The  optimal 
gain  is  found  to  be 


K(k>  = P(kik:  I ) UT(k)  |H(k)  P(k,k  - l)HT<k>  + Rtk)l 
= PtkjkHF  (k)  R 1 tk> 


(2.14a) 

(2.14b) 


Note  that  the  bracketed  term  to  be  inverted  is  of  dimension  <m  X m)  which  is  usually 
smaller  than  the  other  matrices  with  dimension  (n  X n).  The  state  can  then  Ik-  updated, 
using  Equation  (2.8).  ;mu  th>*  ciror  covariance  is  updated  as 


P(k  k)  = (I  - K(k)  lick |)  P(k  k I) 


(2.15) 


We  find  that  the  error  covariance  is  reduced  by  the  update  step,  usually  more  than  it  was 
increased  in  the  extrapola.ion  step,  causing  a net  reduction  in  the  error  covariance  over  a 
complete  -.-ye'e  of  the  filter  (unless  the  filter  has  reached  steady-state  operation).  If  there 
is  no  process  noise  :n  the  system  <Q  = 0).  the  error  covariance  terms  eventually  approach 
zero  (for  the  models  considered  in  this  report).  The  error  covariance  significantly  in- 
fluences the  operation  of  the  filter  throu  ih  (he  gain  matrix,  fn  general,  the  larger  the 
error  covariance,  the  larger  the  gains.  Th  - gain  represents  tin.  relative  weighting  between  the 
old  estimate  and  the  new  data.  Initially  the  error  c'jvariancc  jnd  gain  arc  large,  which 
represents  the  fact  that  we  have  noi  processed  tnucii  data  and  therefore  can  not  place  much 
confidence  in  our  estimate.  At  long  times,  the  error  covariance  and  gains  approach  their 
minimum  values  and  reflect  the  idea  that,  having  processed  a large  quantity  of  data,  we 
place  more  confidence  in  our  estimate.  I*’  there  is  no  process  noise,  eventually  the  error 
covariance  becomes  small  and  the  filter  begins  to  virtually  ignore  any  new  data.  This  con- 
dition can  lead  to  filter  divergence  if  the  dynamic-,  are  - ’re-  suboptiinallv.  This  problem 
will  be  discussed  in  a later  section.  A filter  with  sma5;  ■ • ; covariance  and  srnai!  gains  will 
often  be  referred  to  as  having  a long  "memory”  (or  "narrow  bandwidth”  for  there  who 
consider  filters  from  the  frequency  domain).  Con?.-.-,  ely.  a filter  with  large  error  con- 
variance  and  gains  will  be  said  to  have  a short  memory  or  wide  bandwidth. 

Realization  of  the  optimal  Kalman  filter,  we  have  found,  requires  exact  knowledge  of 
system  dynamics,  measurement  statistics,  and  process  noise  statistics.  In  general,  the  informa- 
tion is  never  available  exactly.  Eoi  example,  the  requirement  for  measurement  statistics 
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poses  some  interesting  problems.  It  is  often  difficult,  even  under  controlled  laboratory  con- 
ditions. to  obtain  good  error  statistics  for  tracking  instruments.  It  is  frequently  even  more 
difficult  to  apply  these  statistics  with  confidence  when  the  sensor  is  employed  in  an  un- 
controlled environment  such  as  a ship  at  sea  operating  against  real  targets.  We  therefore  can 
expect  the  values  of  the  statistics  assumed  for  the  filter  implementation  to  be  slightly  mis- 
matched with  the  real  world.  Such  a mismatch  results  in  a suboptimal  filter  implementation. 
In  order  to  assess  the  effect  of  such  a mismatch,  we  resort  to  a technique  known  as  a 
“sensitivity  analysis.”  whereby  we  calculate  the  actual  covariance  of  the  suboptimal  filter 
implementation  and  compare  it  with  the  optimal  “real-world”  value  possible.  If  a sub- 
optimal filter  of  the  same  structure  as  the  optimal  but  with  unmatched  parameters  (leading 
to  suboptimal  gains)  is  utilized,  the  actual  error  covariance  Pact  can  be  calculated  with  the 
following  equations. 


PACI(k|k-  1)  = 0A(T(k.  k-  l)PA(T<k-  Ik-  1)  0ACT(k,  k - I)  + QACT(k- 

P/U-,  (k|k)  = |l-  K(k)  H A ('T ( k) | PA(  T(klk  - 1 ) [I  - K(k)  HACT(k)] T 
+ K(k)  RACT(k)  KT(k) 


1) 

(2.16) 

(2.17) 


The  subscript  ACT  means  that  the  matrix  is  evaluated  with  the  true  parameters.  The  gain 
K(k)  is  assumed  to  be  calculated  with  assumed  design  values  of  the  parameters.  An  impor- 
tant point  is  that  Kquution  (2.15)  gives  the  actual  error  covariance  only  when  the  optimal 
gain  is  chosen.  Strictly  speaking,  liquations  (2. 16)  and  (2. 17)  can  be  used  only  to  represent 
the  actual  error  covariance  when  a suboptimal  gain  is  chosen  but  the  implementation  of  the 
dynamics  and  measurements  are  assumed  correct.  In  fact,  Gelb  (1974)  (pages  254  and  271) 
argues  that,  under  certain  conditions  present  in  this  application,  these  actual  error  covariance 
equations  are  correct  even  with  incorrect  transition  and  observation  matrices.  In  any  case, 
these  equations,  often  referred  to  as  “simple  sensitivity  equations,”  are  the  equations  used 
for  all  sensitivity  work  in  this  report  and  the  more  complex  “model-reference”  equations 
are  not  used.  Presentation  of  sensitivity  results  has  been  found  to  be  most  recognizable  if 
the  actual  parameters  are  held  fixed  and  the  design  values  allowed  to  vary. 


Since  the  Kalman  filter  equations  are  recursive,  it  is  not  necessary  to  store  a large 
quantity  of  data.  This  factor,  along  with  the  basic  computational  simplicity  of  the  equa- 
tions. allows  easy  implementation  on  a digital  computer.  Notice  that,  since  the  gains  are 
not  a function  of  the  actual  measurements  (for  a truly  linear  system),  it  is  sometimes  pos- 
sible to  precompute  the  gains  and  store  them  as  opposed  to  calculating  them  in  real  time. 
This  precomputation  can  be  performed  only  when  values  of  0(k,  k - 1 ),  R(k),  Q(k),  and 
H(k)  can  be  predetermined  for  all  k.  With  the  initial  value  of  P specified,  all  future  values 
of  P(k)  and  K(k)  can  then  be  calculated.  It  will  be  shown  that,  for  the  Kalman  filter  appli- 
cation in  this  report,  it  is  not  possible  to  predetermine  the  above  matrices  since  they  are 
nonlinear  functions  of  the  actual  target  motion  in  real  time.  The  fact  that  the  Kalman 
filter  must  propagate  the  error  covariance  matrix  as  it  runs  is  an  ambiguous  requirement 
since  the  filter  derives  both  its  power  and  its  principal  computational  burden  from  these 
matrix  equations  for  P. 

The  basic  Kalman  filter  algorithm,  for  the  situation  described  in  this  section  (i.e., 
linear  equations  and  uncorrelated  noise),  is  summarized  in  Table  2. 1 . A FORTRAN  IV 
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Table  2.1.  Kalman  Filter  Algorithm 


Model  x(k  + l)  » p(k  + l.k)x(k)  + w(k)  (I) 

Observations  z(k)  = H(k)x(k)  + v(k)  (2) 

Statistics  E[w(k)|  = E[v(k)|  = Efvfj)  wr(k)l*  = 0 (3) 

E(w(j)wT(k)j  = <Xk)6jk  (4) 

E|v(j)yT(k)|  = R(k)  6jk  (5) 

State  Extrapolation  xtklk-  I)  - #k.  k~  Dxjfk-  lik-  1)  (6) 

Covariance  Extrapolation 

P(kik-  1)  = <Hk.k-  KP(k-  lik-  D^tk.k-  I)  + Q(k-  I)  (7) 

Gain  K(k)  = P(k  k-  l)  HT(k»  {IKUPtk  k - liMr(k)  + Rfk))'1  (8) 

State  Update  xtklk)  = itklk  - l)  + K'k)lz(k)-  H(k)x(kik-  l)J  (‘J) 

Covariance  Update 

Pt'kik)  = |l-  K(k)Htk)}  Ptk.k-  h <IOa> 

= Ptklk-  1)  - K(k)  |H(k) Ptklk - l)HTtk)  + R<k>)  K*(k>  (10b) 


version  of  this  algorithm,  called  KALMAN,  can  be  found  in  Appendix  A.  This  subroutine 
c/as  utilized  for  many  cf  the  simulations  required  for  this  study. 

III.  THE  DERIVATIVE  POLYNOMIAL  TARGET  MODEL 

In  this  section,  the  concepts  and  equations  of  a one-dimensional,  derivative  poly- 
nomial filter  are  introduced.  The  target  model  that  is  ultimately  recommended  (in  Section 
IV)  is  not  exactly  of  a polynomial  form  but  is  quite  similar  in  implementation.  As  such,  the 
parametric  study  of  the  polynomial  filter  in  this  section  is  applicable  ir.  many  ways  to  the 
nonpolynomial  target  model  selected  in  the  next  section  and.  indeed,  is  partly  responsible 
for  that  selection.  It  should  also  be  mentioned  that  the  covariance  convergence  study,  re- 
ported in  this  section,  was  conducted  under  the  assumption  of  optimal  conditions.  By  this, 
we  mean  the  following:  the  actual  target  mod-1  is  always  assumed  to  be  a polynomial  match- 
ing the  filter  model:  the  measurement  errors  are  truly  uncorrdated  and  Gaussian:  and  the 
actual  measurement  error  variances  match  the  assumed  (filter)  measurement  error  variances. 
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Before  continuing,  a very  brief  discussion  of  the  polynomial  target  model  is  in  order. 

In  a general  tracking  filter  application,  where  the  particular  strategy  being  exercised  by 
the  target  is  unknown,  the  form  of  the  state  vector  and  its  propagation  characteristics  3re 
assumed  and  may  or  may  not  represent  the  true  target  motion  over  long  periods  of  time. 

We  do  demand,  however,  the  assumed  target  model  to  provide  a reasonably  accurate  approxi 
mation  of  target  motion  for  short  periods  of  time.  The  choice  of  a polynomial,  of  course, 
is  then  quite  logical  since  the  ability  of  a polynomial  to  approximate  a process-at  least 
locally  is  well  known.  Indeed,  the  famous  Weierstrass  approximation  theorem  tells  us  that 
polynomials  can  approximate  an>  continuous*  function  over  a finite  interval  to  any  degree 
of  approximation  desired.  We  are  not  really  interested  in  considering  very-high-order  poly- 
nomials. however.  It  must  be  emphj sized  that  the  gunfire  control  problem  is  not  primarily 
a fitting  or  smoothing  problem  but  a prediction  problem.  A higher-order  (and  thus  more 
accurate)  po./nomia!  fit  of  observed  past  data  may  not  help  us.  and  may  even  hinder  us.  to 
more  accurately  predict  future  target  po-ition  for  long  extrapolation  times.  Also,  higher 
order  filters  are  computationally  more  burdensome  and.  dcpcndii  g on  >he  order  of  the 
highest  derivative,  may  not  settle  fast  enough  for  the  gunfire  control  problem  (as  will  be 
diown).  For  these  reasons,  we  consider  polynomials  with  nonvanishing  derivatives  only 
through  third-order  (jerque).  The  important  advantages  of  such  low-order  polynomials  are 
computational  efficiency  and  limited  necessity  for  n priori  target  maneuver  strategies.  For 
a complete  treatment  of  various  polynomial  filters,  smoothers  and  predictors,  the  reader  is 
referred  to  Morrison  (1969). 


A.  KQUATIONS 

Derivative  polynomial  models  are  described  by  (he  differential  equation 

d"*1  x(t)/dt»*‘  = 0 (3.1) 

where  dn  x/dtn  is  the  highest  nonvanishing  derivative  and  n :s  referred  to  as  the  ‘order”  of 
the  model.  The  state  vector  b defined 

x(t)  = fx(t)  dx(t)/dt  d-x(t)/dt2 

...  dnx(t)/dtn  JT  (3.2) 

= | Xj (t)  Xitt)  Xj(t)...  Xn*|(t)JT 


* All  unccl*  *tih  mja  mu't.  of  coane,  be  con(im*<«n  k>  portion  and  wioafy  Acotlcrjlkw  ma>!  jlwayt  be  Out:  bat  may 
Nr  efTei, lively  diKonimuo  iv  Tbiv  implies  yctquc  may  sol  always  cxisl. 
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Notice  that  the  state  vector  has  dimension  (n  + 1 ).  In  the  notation  of  continuous  linear 
systems. 


dx(t)/dt  = Fft)x(t) 


We  can  write  the  elements  of  the  matrix  F(t)  as 


fij(t)  = &i.j-  i I < i-j  < (n+  !» 
where  5,j  | is  the  Kroneckcr  delta.  The  transition  matrix  is  therefore 


where 


0,.(k  + 1.  k)  = (ADi^/tj-i)!  1 < i.j  < (n+  I) 


At  = tlk  + I)  - tck»  = I/C 


l/iji).  = 0 for  j < i 


For  example,  for  n - 2.  we  find 


otk  + I.  k>  = 0 = 0 ! 


0 0 


l At  At*’/2 


For  purposes  of  this  report,  the  data  rate  k tiler t/t  will  always  be  assumed  to  be  constant  so 
that  the  transition  matrix  (a  function  of  At  only)  is  time-invariant.  For  the  present,  we  will 
assume  process  noise  is  zero. 

In  our  application,  we  will  have  available  measurements  of  position  oniy.  so  that  the 
measurement  vector  and  covariance  reduce  to  scalars.  That  is. 

/(k)  = | /<  k ) j <3.9) 

R(k)  = Jo’tkti  <3.10) 


where  o-  is  the  measurement  error  variance.  The  elements  of  the  (ni  X tn  ♦ I t)  observation 
matrix  are  simply 


lijjdO  = 5j|  5ji  i = m = I 

I = j *S  (n  + 1) 


C5.ll) 


For  n = 2,  we  find  the  row  vector 


H(k)  =!!=(!  0 0] 


The  residual  also  reduces  to  a scalar  quantity 


Kklk-  I)  = zCk)  - x i ( k |k  - I) 


(3.12) 


(3.13) 


Hie  simplicity  of  the  measurement  matrices  result  in  the  gain  calculations  being  effected  as 
rather  simple  algebraic  functions  of  the  extrapolated  error  covariance.  K.g.,  the  gain  is  the 
column  vector  with  elements 

Kj(k)  = 1’ijfklk-  D/JP, , (k|k  - l)  + o2(k)|  j = 1,  n+l  (3.14) 


I).  INITIALIZATION 

In  order  to  specify  completely  a set  of  equations  representing  a derivative  polynomial 
Kalman  filter,  we  have  to  consider  a method  of  initialization.  Since  the  Kalman  filter  equa- 
tions are  recursive,  values  of  both  the  state  and  associated  error  covariance  from  the  previous 
cycle  are  required.  For  most  linear  filters,  the  method  of  initialization  is  not  particularly 
significant  because  the  effects  tend  to  vanish  at  large  times.  Indeed,  that  is  the  case  for  this 
problem.  However,  it  is  desirable  to  limit  the  excursions  of  the  predictor  while  the  filter  is 
in  the  unconverged  condition.  One  basically  wants  to  form  the  best  state  estimate,  based 
on  any  information  that  is  available,  that  minimizes  the  associated  initial  error  covariance. 

In  the  absence  of  Any  a priori  target  information,  the  only  data  available  to  initialize  the 
filter  will  be  the  measurements.  In  order  to  initialize  the  derivative  elements  of  the  state 
vector,  it  is  necessary  to  have  available,  at  a minimum,  enough  measurements  to  define  the 
required  derivatives  by  differing  techniques.  This  method  approximates  derivatives  (known 
as  ‘‘pseudo-measurements”)  by  finite  differences  of  current  and  prior  position  measurements. 
For  example,  let  us  consider  the  case  n = 2 and  use  the  well-known  backward  difference 
formulae  (neglecting  higher-order  terms)  where  we  approximate 


x,  (010)  = z(0) 

x2(0|0)  = |z(0)~  z(-l)|/At 

x3(0|0)  = | z(0)  ~ 2z.(-  I ) + z(-2)|/At2 


(3.15a) 

(3.15b) 

(3.15c) 
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We  can  now  evaluate  Pt'O/'O)  which  we  will  write  in  the  normalized  form 


Pfkik)  - o,(k) p(k) o, (k) 


(3.16) 


where  o,(klk)is  the  diaganol  matrix  of  standard  deviations  of  the  state  vector  error  Ojfkik) 
and  p(k)  is  the  matrix  of  normalized  correlation  coefficients.  Assuming  the  measurement 
error  variance  is  not  rapidly  changing  and  neglecting  finite-difference  approximation  error, 
we  find 


3,(010)  = o(0)  0 </l!A( 


0 I 


(3.57a) 


p(0)  = ! I.V- 


lfy/b 

x/?'- 


_1  - 1 | 


Unfortunately,  for  high-data-rate  systems,  the  error  variances  of  the  derivatives,  using  the 
finite-difference  method  of  state  estimation,  can  become  unreasonably  large.  For  example, 
from  Hquation  (3. 1 7a).  we  find  that  for  o(0)  = 5 yards  and  Z)t  = 1/16  second,  the  standard 
deviation  of  the  acceleration  error  is  approximately  2lf2  G's!  This  situation  is  clearly  un- 
desirable since  we  could  guess  zero  for  acceleration  and  never  realize  an  error  larger  than  a 
few  G’s.  A similar  argument  can  be  made  concerning  the  velocity  initialization  but  the 
numbers  are  not  as  dramatic.  Therefore,  we  will  choose  a single-pass  (using  one  piece  of 
data)  initialization  for  the  second  order  filter  bv  taking 

x,  <0!  0)  = z(0>  (3.18a 

Xi(0l0)  = 0 (3.18b 

x t (01 0)  = 0 (3.18c 


Thc  values  for  the  derivative  state  error  standard  deviations  would  then  reflect  our  estimates 
of  the  possible  errors  That  liquation  (3. IS)  would  produce.  This  is  effected  3s 
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0 


0 


where 


a'-iO) 


PfOlO) 


0 01(0)  0 

0 0 o|{0) 


oi(0)  = !•:  |dx(0)/dt  j - 
05(G)  = !■:  I d- x*0>/dt-  ) - 


o.m 


(3.20a) 


(3.20b » 


area  priori  estimates  and  would  depend  upon  the  type  of  target.  For  example,  it  has  been, 
found  that  o»( 0)  = 300  yards/second  (approximately  Mach  f ) and  oyiO)  = 10  yards  second- 
t approximate!}  one  C)  quite  adequately  span  the  spectrum  of  possible  air  targets  i:i  the  cur- 
rent target  scenario.  Singer  ( !‘>70)  has  suggested  a method  of  estimating  o ;(0l  as 


o^lO)  = — ~r— - }l  + 4 • Probability  i I x 1 1 = Amax ) - Probability  (xs  = 0)j  (3.2!) 

where  AmjX  is  the  maximum  acceleration  capability  of  the  target  which,  presumably,  one 
might  estimate  for  a particular  class  of  targets.  The  values  of  oWO)  and  03(D)  for  surface 
targets  would  undoubtedly  be  much  smaller.  For  example,  we  might  choose  a2 (0)  = 

14  yards/second  (about  25  knots)  and  o»(0)  - 2 yard  second-  tabout  0.2  (i).  Upon  com- 
paring error-convariance  results  of  the  one*p:iss  initialisation  with  those  of  the  finite  dif- 
ference initialization,  it  was  found  that,  after  approximately  10  seconds,  there  was  no 
appreciable  difference.  On  the  other  hand,  the  one-pass  initialization  yields  considerable- 
lower  error  during  this  initial  period.  For  the  case  of  n = 3.  when  a value  for  04(G)  is  re- 
quired. we  must  estimate  the  maximum  rate  of  change  of  acceleration.  For  a high- 
performance  air  target,  such  as  a modem  fighter  aircraft,  the  author  learned  (from  dis- 
cussions with  pilots)  that  the  maximum  sicadv-sian-  turning  rate  over  the  entire  speed/ 
maneuverability  spectrum  is  approximately  20 degrees/second  = 0.3-1')  radians/second.. 
A value  of 

o4(0)  = u }„  * o}(Qi  (3.22) 

was  therefore  chosen  for  initialization  of  the  third  order  niters. 

In  conclusion,  the  one-pass  initialization  is  very  simple  to  implement  as  it  requires  no 
additional  measurements  or  logic.  It  was  found,  both  from  error,  covariance  results  and 
Monte  Carlo  simulation,  to  yield  improved  estimates  during  the  initial  transient  phase,  it  is 
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used  in  ail  remaining  work  in  this  report  and  was  subsequently  implemented  in  the  digital 
MARK  68  GFCS. 

C.  CONVERGENCE  PROPERTIES 

Early  in  the  course  of  this  investigation,  it  was  discovered  that  fundamental  information 
on  factors  that  directly  affect  filter  convergence  and  settling  time -and  therefore  filter 
performance - was  not  available.  Qualitative  relationships  were  usually  known,  or  at  least 
suspected,  between  convergence  rates  and  such  variables  as  polynomial  order,  data  rate, 
initialization  methods,  etc.  It  was  felt,  however,  that  quantitative  measures  of  such  effects 
were  necessary  if  intelligent  decisions  and  tradeoffs  v/ere  to  be  made  concerning  the  basic 
operating  conditions  of  the  filter.  As  mentioned  previously,  these  studies  v/ere  conducted 
under  optimum  operating  conditions  and.  as  such,  are  intended  to  reflect  relative  perform- 
ance  of  one  filter/parameter  set  with  respect  to  another.  • 

The  two  significant  parameters  that  describe  a polynomial  filter  are  the  polynomial 
order  n and  the  data  (or  cycle)  rate 

dk-dt  = k = 1 M 

In  order  to  establish  the  relationships  between  filter  performance  and  these  parameters,  the 
filter  matrices  of  the  preceding  paragraphs  were  inserted  into  the  general  Kalman  filter 
subprogram  in  Appendix  A which  was  exercised  for  various  values  of  these  parameters. 
Specifically,  the  values  n = 1.2.  and  3.  corresponding  to  constant-velocity  (first-order), 
constant-acceleration  (second-order)  and  constant-jerque  (third-order),  were  selected  Data 
rates  of  k = 2.  4.  X.  1 6.  and  32  Hertz  were  considered  to  span  the  range  of  reasonable  pos- 
sibilities. Assuming  a t>  pica!  constant  measurement  error  standard  deviation  and  zero- 
process  noise  Q.  the  error  covariance  and  gains  can  be  calculated  independently  and  no  target 
measurements  need  be  simulated.  All  conditions  were  assumed  to  be  optimal  so  that  filter 
error  >.o\aiiance  would  truly  represent  filter  convergence.  In  keeping  with  the  idea,  pre- 
viously expressed,  that  predicted  position  is  the  important  quantity,  the  quantity  actually 
used  by  a fire  control  system,  the  standard  deviation  of  predicted  position  error  o\  (t  + 
tpit)  was  chosen  as  the  ultimate  measure  of  relative  covergenee.  Piedieted  position,  of 
course,  can  not  settle  until  all  elements  of  the  state  vector  have  settled  so  that  0|  (t  + 
tp  t)  is  an  effective  '’norm'*  of  the  current  error  covariance.  In  fact,  it  can  be  shown  that 
any  cost  function  of  the  form 

J = !•:  |e'  A e| 

can  be  uniquely  minimized  with  respect  to  K only  if  A is  positive  semi-definite  (see  Cleib. 

I ‘>75  for  exampie).  Of  course,  the  trace  of  P is  only  a special  case  when  A = I.  One  could 
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also  choose  the  positive  definite  matrix  A = 0T(tp)  0(tp)  as  we  have  and  minimize  the  pre- 
diction error.  For  the  case  of  n = 2,  we  are  interested  in 


o,(t  + tpit)  = (P|, (tit)  + 2tp  Pi 2 (t ! t)  + tp  P22(tlt)  + tJ  Pi3(tlt) 
+ tj>P2.  ■t)  + tJP33(tit)/41'/2 


(3.23) 


A prediction  time  of  tp  = 10  seconds  is  used  as  a standard  of  comparison  for  all  the  one- 
dimensional results  in  this  study.  This  corresponds  to  a nominal  projectile  time-of-flight  for 
a Navy  five-inch  gun  system  engaging  air  targets.  Notice  also  in  Equation  (3.23)  the  factors 
that  multiply  the  rate  errors.  It  can  clearly  be  seen  that  velocity  and  acceleration  estimation 
errors  are  the  principal  cause  of  prediction  noise. 

Before  discussing  the  convergence  results,  it  might  be  helpful  to  consider  a transforma- 
tion of  the  polynomial  filter,  suggested  by  Morrison  (1969),  that  is  useful  in  understanding 
and  interpreting  these  results.  A transformation  of  the  state  vector  can  be  used  to  absorb 
the  At  dependence  in  the  error  covariance  equations.  This  is  done  by  defining  an  alternate 
state  vector 


x*  (t)  - [x  (At)dx/dt  (At2 /2) d2 x/dt-  .... 

...  <Atn/n!)d"x/dtnjT 


(3.23) 


The  equivalent  transition  matrix  is 

0i*  = (j-  DVf(i-  D!(j-  i)!|  (3.24) 

which,  we  find,  is  a matrix  of  constants  independent  of  At.  For  the  example,  n = 2.  we  find 

■'  1 <1 


0’  = 


0 


0 0 


It  the  At  dependence  is  removed  from  the  initial  error  covariance,  the  error  covariance  a.-d 
gains  are  then  a function  only  of  the  data  point  k and  the  order  n.  In  subsequent  sections, 
we  will  show  that  the  frequency  content  of  target  motion  is  ratbei  low  so  that  a high  sam- 
pling rate  is  not  actually  necessary  to  satisfy  the  Shannon  sampling  theorem.  In  fact,  the 
principal  advantage  of  a high  data  rate  is  to  reduce  the  effective  measurement  error  level.  In 
the  section  on  prefiltering,  it  will  be  shown  that  data  compression  techniques  can  be  used 
to  achieve  effective  high  data  rates  without  actually  cycling  the  Kalman  filter  at  such  rates. 


18 


fn  order  to  consider  the  real  convergence  criterion,  the  extrapolated  position  error 
standard  deviation,  it  is  more  straightforward  to  use  the  original  unsealed  equations  (without 
Morrison  s transformation).  In  Figure  3.1 , we  have  plotted  the  time  history  of  the  standard 
deviation  of  the  predicted  position  error  (tp  = 10  seconds),  normalized  by  a (the  measure- 
ment error  level),  for  each  filter  order  and  data  rate  pair  (n.  &)  under  consideration.  Such  a 
graph  can  be  used  to  determine  optimal  settling  time  for  any  particular  situation.  For 
example,  suppose  we  decide  the  filter  is  “settled"  enough  to  fire  the  gun  when  the  predicted 
target  position  is  accurate  to  within  25  yards.  If  we  can  track  the  target  with  a a = 5 yards, 
then  using  the  (2.  16)  filter,  we  could  start  firing  8 seconds  after  the  initiation  of  target 
track.  Basically,  we  see  from  this  graph  that  the  lower  the  order  and  the  higher  the  data 
rate,  the  faster  convergence  will  occur.  We  see  that  the  third-order  filter  converges  painfully 
slowly  at  these  data  rates  and  is  probably  not  suitable  for  gunfire  control  purposes  due  tn 
this  factor  and  the  heavy  real-time  computational  burden  of  this  filter  at  high  data  rates 
The  constant-velocity  filter,  on  the  other  hand,  converges  very  rapidly,  predicting  better 
than  the  tracking  measurement  error  in  less  than  10  seconds,  with  a data  rate  of  only  4 Hertz. 


Of  course,  settling  time  of  these  optimal  filters  is  not  the  only  consideration  in  making 
a choice.  We  know  that  targets  do  accelerate  and  that  the  constant-velocity  filter,  without 
an  estimate  of  acceleration,  will  do  quite  poorly  in  predicting  future  target  position.  The 
error  in  that  situation  is  of  a bias  type  rather  than  the  uncorrelated  error  considered  here. 

In  Figure  3.2.  we  have  plotted  for  reference  purposes  the  complete  error  convergence  history 
of  all  the  state  elements  (position,  velocity,  acceleration,  and  1 0-second  predicted  position) 
for  the  constant  acceleration  (n  = 2)  filter.  With  the  dashed  lines,  we  have  also  plotted  the 
corresponding  actual  error  that  would  occur  if  we  attempted  to  filter  the  “real-world" 
constant  acceleration  model  by  assuming  a constant  velocity  (n  = 1 ) filter  model.  I.arge 
bias  errors  in  ail  state  elements  are  observed  to  occur  in  the  actual  error  elements  when  using 
an  n = 1 filter  in  an  n = 2 environment.  The  prediction  error,  which  initially  starts  to  get 
smaller,  quickly  reverses  itself  and  satisfactory  prediction  error  is  never  achieved.  The  n - 1 
filter  i<  really  not  a candidate  if  we  are  going  to  attempt  to  deal  with  manuevering  targets. 
That  leaves  the  n = 2 filler  which  we  would  choose  to  implement  if  a polynomial  target 
model  were  selected.  We  would  also  choose  the  k = 16  Hertz  (or  faster)  since  it  converges 
fastest.  Convergence  time  is  not  only  important  initially,  but.  as  we  shall  see  in  a late,  sec- 
tion on  adaptation,  it  may  be  necessary  for  the  filter  to  reconverge  when  the  adaptive 
tea  turns  in  the  filter  decide  it  necessary.  Filter  responsiveness  in  such  a situation  dictates 
greater  significance  to  the  problem  of  realizing  a short  settling  time. 


In  conclusion,  then,  we  would  choose  the  second-order  filter  over  the  third-order 
(ruled  out  due  to  settling  time)  and  over  the  first-order  (because  we  must  track  an  accelerat- 
ing target  in  an  effective  manner).  We  also  conclude  that,  given  uncorrclatcd  measurement 
noise  (satisfying  the  requirement  for  the  derivation  of  this  filter),  we  want  to  process  data 
at  the  highest  rate  possible.  The  finai  determination  of  this  rate  will  be  based  on  available 
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Figure  3.1.  Derivative  Polynomial  Fillers 
with  Uncorrclated  One-Pass 
Initialization  Normalized 
Standard  Deviation  of  Pre- 
dicted (10  sec)  Position  Error 


Tim* 


Figure  3.2.  Optimal  Convergence  of 
n = 2 Filter  and  Actual 
Error  Using  n = 1 Filter 


computation  time  in  the  real-time  filter  implementation  and  on  the  ability  to  determine  and 
remove  the  actual  correlated  portio.:  of  the  measurement  no»sc.  This  la:i  aspect  of  the 
problem  will  be  discussed  in  a later  section.  We  will  pircced  on  the  basis  of  k = 16  Hertz 
unless  otherwise  noted.  For  reference  pumoses.  a summary  of  the  filter  equations  for 
n = 2 can  be  found  in  Table  3. 1 . 

Another  possibility,  which  the  author  would  like  to  exrlort  in  the  future,  is  that  of 
polynomial  fitter  cascading.  In  this  approach,  one  could  conceivably  obtain  *bc  fast  con- 
vergence properties  of  the  low-order  filter  and  the  modeling  improvement  of  the  high-order 
filters.  The  procedure  might  be  implemented  in  the  following  manner.  A constant  velocity 
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filter  (n  = 1).  with  limited  memory  to  prevent  divergence  (t  be  discussed  lute:-).  would  first 
be  implements.  rl.  The  output  state  v ctorof  this  filter  would  then  provide  the  measure- 
ments 3tisi  measurement  error  cr-  .ianee  for  a constant-acceleration  tr>  ~ 2)  fitter  false  with 
specified  memory).  It  might  even  hr  p rssthrV  to  us?  another  cascade  and  feed  a jerque 
in  = 3>  filter.  A logic  section  would  then  operate  on  the  residuals  of  each  filter  and  pick  the 
one  tbft  is  operating  the  best.  The  possibilities  are  certainly  iniriguhts  ~v.i  should  he  ex- 
plored if  possible. 
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Tabic  3.1.  Second-Order  Filtei  bqustions:  A Summar** -(Continued) 


S 

K,(k)  = Pn(klk-  t)/|P||(klk-  !)  + o2(k)j  (10) 

K2(k)  = P,2(kik-  l)/(P,,(k!k-  !)  + o2(k)l  (11) 

Ki(k)  = P,  j(k!k-  D/iPjjtkik  - ?)  + o2(k))  (12) 

KHSIDUAI. 

rlkik  - i)  *'  4tK,  - i (kik-  I) 

CPDATI: 

x«  ikiki  = 'i  -tkik  - » » + K| <k) »--(kik -•  I ) 

\*;k-'k)  - ( k i k - !)  + Kiiklf-tkik - 1) 

\jtk.kJ  = \j(klk-  l)  + Kitk-i-lkik  - i 

P,,«k.k>  * II  - K i <k)i  P, | (klk  - ! ) 

P,:«k!k>  * {1  - kj (k)j  P,:(klk~  t) 

Pi?(kik)  = {I  - K|<k)i  Pj j(kik-  I) 
p»;iklk»  = Pt>(klk-  1 )~  k2(k) P| 2< k i k - 1) 

P>?(kik)  * P>.n«k  tk  - I)  - k;<k)Pi3(kik-  1) 

P3 .; ( k i k > *-  Ptjiklk - »)  - K3(k)Pn(kik-  1) 

INITIALIZATION 
x | (OK))  = /.(0» 
x?(0io>  = \»(Oio>  = o 


Tabic  3.1.  Second-Order  Filter  Equations:  A Summary~(Continued) 


INITIALlZATION-(Continued) 

P,,(0i0)  = o2(0) 

(25) 

P,  2(0 10)  = 0. 

(26) 

P{3«0l0)  = 0. 

(27) 

P22(OlO)  = 5?(0) 

(28) 

P23(OIO)  = 0. 

(29) 

P33(010>  = oj(0) 

(30) 

PREDICTION 

x,  (t  + tpil)  = x,  (tit)  + x 2 (tit)  tp  + xj(tltp2/2 

(31) 

RECOMMENDED  VALUES  OF  PARAMETERS 

k = 16  Hertz 

(32) 

FOR  AIR  TARGETS 

o2*0)  = 300  yards/second 

(33) 

O3(0»  - 10  yards/second2 

(34) 

FOk  SURFACE  (' . IP)  TARGETS 

(35) 


O2(0)  = 14  yards/second 
oj(0)  = 2 yard/sccond2 


(36) 


IV.  Till:  RANDOM-ACCELERATION  TARGET  MODEL 


The  derivative  polynomial  target  model  of  Section  III.  while  possessing  several  advan- 
tages. undoubtedly  suffers  as  far  as  being  a realistic  representation  of  the  scenario  of  actual 
targets  over  a period  of  more  than  a few  seconds.  Actual  targets  would  rarely  traverse  a 
simple  predictable  path  while  in  the  range  of  engagement  by  our  GFCS.  Instead,  such 
targets  would  be  mamicrerin g.  These  maneuvers  could  take  the  form  of  a deterministic  path 
if  the  target  has  a particular  goal  in  mind,  such  as  effecting  a collision  with  the  ship  (in  the 
case  of  an  anti-ship  missile)  or  the  release  of  ordnance  (in  the  case  of  a bomber).  The 
scenario  would  also  have  to  include  nondetenninistic  (or  random)  maneuvers  in  such  situa- 
tions w-her.  the  target  is  maneuvering  merely  to  evade  interception  by  the  ship's  weapons  or 
when  irregular  winds  and  atmospheric  turbulence  would  act  on  the  target.  Anti-ship  missiles 
would  be  subject  to  the  latter  effects  and  also  may  respond  to  the  wander  of  its  own  radar 
around  the  target  ship.  In  any  ease,  the  probability  of  a target  following  a constant  accelera- 
tion (or  any  polynomial ) course  for  any  length  of  time  is.  unfortunately,  not  what  we  might 
prefer.  With  this  idea  in  mind,  perhaps  the  problem  can  be  approached  from  a different 
point  of  view.  Rather  than  attempt  to  model  the  target  with  a deterministic  trajectory  such 
as  a polynomial,  let  us  consider  a case  where  we  acknowledge  the  fact  that  the  target  is  very 
likely  to  be  maneuvering  in  some  unknown  manner  and  where  we  assume  that  each  dimen- 
sion of  such  trajectories  appear  (to  the  FCS>  to  be  a random  variable.  The  following  model 
was  developed  by  Singer  1 1970)  and  has  already  been  selected  for  implementation  in  several 
other  fire  control  and  tracking  systems. 


The  target  acceleration  is  assumed  to  be  an  autocorrelated  (or  serially  correlated) 
random  variable  with  known  statistics.  Specifically.  Singer  chose  to  model  the  acceleration 
as  a hrst-order  Markov  process  with  the  differential  equation 


/d-x(t)\ 

JL 

/d:x(t>\ 

\ dt:  j 

l «>  ) 

driven  by  a white  noise  input  of  variance 


o£ (T)  = 


The  acceleration  sequence  generated  is  a form  of  process  known  as  a random  walk.  An 
example  of  such  a process  generated  on  the  computer  with  a correlated  random  number 
generator  CORNUM  (See  Appendix  A)  is  shown  in  Figure  4. 1 . Notice  that  the  correspond- 
ing position  (double  integral  of  x)  is  also  plotted.  The  lutocorrelation  function  for  the 
ac«.eler::!ioi.  is  the  exponential 
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Figure  4.!.  Example:  Random  Acceleration  Model 
(rm  - 20  seconds.  om  = 2 yards/second2 ) 

n^„(T>  = hi  xi  i)  x(t  + Tii 

= o*,  exp  ( - iTI/rm) 


We  therefore  require  only  two  parameters  (on.  and  rm  or  tom  ) to  characterize  the  target 
maneuverability.  The  parameter  om  can  be  (loosely)  thought  of  as  the  root-mean-square 
level  of  target  acceleration  and  rm  as  the  characteristic  maneuver  time  (approximately  0.8 
times  the  mean  time  between  acceptation  zeros)  or  the  inverse  of  maneuver  frequency  com 
Notice  that  since  om  is  essentially  the  rms  level  of  acceleration  for  one  dimension,  the 
actual  peak  total  (three-dimensional)  acceleration  might  be  several  (perhaps  three  to  four) 
times  larger  than  om  . It  is  felt  that  this  model  is  much  better  suited  to  our  purposes  than 
the  polynomial  models  because  it  directly  relates  the  filter  parameters  to  the  kinematics  of 
the  target  scenario.  We  have  also  found  that  this  model  performs  better  in  both  filtering 
and  prediction  than  either  the  n = I or  n * 2 polynomial  models.  It  should  be  emphasized 
that  we  do  not  mean  to  imply  that  one  cannot  find  particular  trajectories  that  might  poorly 
rnatch  this  model.  We  do  expect  actual  trajectory  to  be  "contained’’  statistically  in  the 
random  acceleration  model  and  thus  be  a "reasonable”  sample  from  our  assumeu  population. 
The  author  feels  that  the  random  acceleration  model,  while  far  from  perfect,  is  the  best 
mode!  yet  developed  for  general  tracking  filter  application.  The  model  can  be  written  in 
state  space  notation  by  defining  the  state  vector 
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where 


X<t)  = |x(l»  dx(t)/dt  d2xtt)/dt2| 


dx(t)/dt  = Ftt)x(t)  + wit) 


F(t)  = F = 0 0 


0 0-  l/rmj 


wit)  =|0  0 wit)}1  14. 

By  simple  integration  of  liquation  (4.6).  we  find  that  the  transition  matrix  can  be  written 


i At  of  At ) 


04k.  k - I)  = 0 = 10  I ptAl) 


0 0 TtAt) 


where 


TtAt)  = expt-  At irm  ) 

0tAt)  = rm  (I  - 7(At)J 
atAt)  = t2  IHAt)  + At/rw  - !| 


The  argument  of  a.tf  -nd  y - At  - will  hereafter  be  suppressed  unless  specified  for  a par- 
ticular pun»ose.  The  resemblance  of  this  model  to  the  second-order  polynomial  model  can 
be  seen  immediately  if  we  consider  the  limiting  case  of  rm»  At.  Ti.'  transition  matrix 
becomes: 


I 


r 


I 


$ — $2  ~ 


0 


i m At  2n 


At 


0 


(4.10) 


i.e,  identical  to  the  second-order  polynomial  model.  On  the  other  hand,  if  rm  vanishes  (or 
rm  « At  ),  we  find 


I 


At  0 


0 


1 0 


(4. i 1 ) 


0 0 0 


This  transition  matrix  corresponds  to  the  first-order  polynomial  model.  We  would  tliereforc 
expect,  finite  values  of  rm  (between  zero  and  infinity)  to  yield  a model  that  exhibits  con- 
vergence properties  between  those  of  the  first-  and  second-order  polynomial  filters.  Indeed, 
this  turns  out  to  be  the  case.  In  Fipire  4.2.  the  extrapolated  position  error  standard  devia- 
tion is  again  plotted  for  various  values  of  rm . The  cases  rm  = 0 and  rm  -*•  <»  corresponds 
exactly  to  the  first-  and  second-order  polynomial  cm's.  All  cases  were  run  with  a data  rate 
of  16  Hertz,  a value  of  o„,  =0  (limiting  case  of  am  small)  and  with  one-pass  initialization. 
We  find  that,  just  as  expected,  as  rm  gets  larger,  the  convergence  properties  move  toward 
the  second-order  polynomial  case. 


The  process  noise  matrix.  Q.  is  developed  in  the  appendix  of  Singer’s  paper.  This 
development  is  rather  lengthy  and  will  not  be  reproduced  here.  The  final  results  can  be 
written 


On  = r*m  1 1 + HAt/T*  ) - 2(At/rm  )'- 

(4.12a) 

+ 2(At/rm  )’/3  - 4(At/rm  ) y - y2 1 
Qt  2 = T}m  ai  ( I - 72  - -7  + -(A(/rw ) y 


- 2(At/rm  ) + (At/rm  )J  | 

On  * t-  o^,  ( I - y2  - 2(At frm  >7)  (4.1 2c) 

022  = ti  oj",  (4y-  3-  y2  »■  2(At/rm )|  (4.l2d) 

023  * rm  0*  (I  +72  - 27)  (4.1 2c) 

O33  = erj,  ( ! - 72  > **.120 
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Figure  4.2.  Random  Acceleration  Filter 

with  Various  Values  of  normal- 
ized Standard  Deviation  of 
Predicted  < 10  see)  Position 
Error 
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(Notice  that  Q is  the  only  place  where  a^  appears  in  the  filter.)  Again  for  the  limiting  case 
of  ry  > At.  this  matrix  can  be  approximated  as 


AH/20  AH/8  At2  lb 


Q =2At  o*/rM  At3/8  At-/6 


(4.13) 


AH/6 


Since  we  will  always  operate  with  parameters  At  and  r*  in  this  latter  range,  we  will  use 
this  approximate  form  of  Q for  our  work.  Simulations  show,  that  for  values  of  the  parameters 
At  and  in  the  ranges  considered  in  this  report,  there  is  no  difference  between  the  use 
of  Equation  (4. 1 3)  instead  of  Equation  (4. 1 2).  We  must  use  the  exact  expression  for  Q in 
the  long-time  extrapolated  covariance  however.  Notice  in  Equation  (4.13)  that  a decrease 
in  I'm  or  an  increase  in  frequency  influences  Q in  the  same  way  as  an  increase  in  ojj, . The 
extrapolated  position  variance  equation  for  this  case  is 

P|»ft  + t?/t)  " P|  | Ct/t)  + 2 tp  P,  »(t/t)  + 2 Op  Pi  j(t/t) 

+ tj  p2»(t/t)  + 2 Op  tp  P23(t/t)  + ot;  P3J(t/t)  (4.14) 


+ Olip 


I 


where 


Op  = rj,  l exp  (-  tp/Tji)  + tp/tM  - ! | <4.1 5) 

and 

Qnp  = ri  °i«  ! 1 + -ftp/** ) - -(»-'**  )3 

+ 2(tp/ry)5/3  - 4(tP/T||  >exp(~  tp/fy  » (4,16) 

- exp  (-  2 tp/T*)! 

The  choice  of  a good  value  of  ry  b a function  primarily  of  the  target  scenario  for  the 
GITS.  The  autocorrelation  was  studied  for  several  typical  profiles  of  older  anti-ship  mis- 
siles. and  a value  of  approximately  ry  = 20  seconds  appeared  to  consistently  be  yielded. 

Singer  ( 1 970)  recommends  ry  = 20  seconds  also  for  manned  maneuvering  targets  exercis- 
ing evasive  maneuvers.  The  same  value  has  also  been  found  independently  by  other  people 
who  have  studied  the  problem.  No  information  cn  maneuvering  surface  targets  has  yet  been 
analyzed.  Other  values  can.  of  course,  be  chosen  for  surface  targets  or  to  reduce  settling  time 
as  may  be  required.  It  should  be  emphasized  that  this  rather  low  observed  maneuver  fre- 
quency. corresponding  to  ry  = 20  seconds,  does  not  necessarily  imply  that  this  is  the  high- 
est frequency  that  a particular  target  might  be  capable  of  sustaining.  On  the  contrary,  most 
air  targets  can  maneuver  much  more  rapidly  if  desired.  Rather,  a low  maneuver  i requeues 
is  probably  typical  of  air  targets  maneuvering  to  achieve  a particular  goal  such  a*  inter- 
cepting ownship  or  exercising  their  own  fire  control  in  order  to  release  ordnanct  at  own-  s 

ship.  We  consider  values  of  ry  anywhere  in  the  range  of  3 to  20  seconds  as  realistic. 

The  selection  of  a value  of  Oy  is  more  difficult  than  that  for  ry . The  acceleration 
autocorrelation  study  of  the  target  scenario  demonstrates  a very  wide  range  of  values  for 
oy . For  example,  some  targets  might  achieve  an  rms  maneuver  level  dose  to  one  G while 
others  are  essentially  nonmaneuvering  (actually  maneuvering  at  a very  low  level  due  to 
atmospheric  effects).  The  effect,  however,  of  o^  on  optimal  filter  performance  is  quite 
dramatic.  For  example,  in  Figure  4.3.  the  steady-state  value  of  the  normalized  10-sccond 
prediction  error  is  plotted  as  a function  of  Oy  for  ry  = 20  seconds.  For  a value  of  oy  = 

10  yards/second-  (at  the  upper  limit  close  to  one  G).  the  normalized  prediction  error  is 

15.96  . clearly  a very  poor  value  for  GFCS  applications  with  conventional  projectiles.  The  J 

reason  such  a large  prediction  error  occurs  is  due  to  the  value  of  the  current  error  covariance 
which  remains  large  due  to  the  process  noise. 

Figure  4.3  also  plots  the  current  acceleration  estimation  error  level  which  we  see  is 
quite  sensitive  to  o y . Since  the  prediction  error  equation  (4.14)  magnifies  the  acceleration 
error,  we  find  the  large  prediction  errors  occurring  for  large  oy  . Obviously,  the  parameter 
oy  has  a significant  effect  on  the  filter  bandwidth  and.  since  it  has  such  a wide  range  of 
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values,  must  be  carefully  chosen  to  achieve  proper  filter  operation.  The  next  section  on 
adaption  will  deal  in  more  detail  with  the  specification  of  Oy  and  ry  • For  reference  pur- 
poses. the  convergence  of  the  normalized  prediction  error  is  plotted  in  Figure  4 .4  for  a 
range  of  values  of  oM  and  vM  . Notice  that  the  steady  state  values  for  smaller  rM  (higher  ) 
are  nut  necessarily  larger  for  prediction  since  smdrer  Ty  suppresses  prediction  noise  at  the 
same  time  it  increases  filler  noise. 

It  is  also  interesting  to  determine  the  sensitivity  to  incorrect  assumed  values  for  the 
random  acceleration  parameter.  To  do  this,  we  first  determine  a set  of  parameters  Oy  f and 
rmS  tor  the  filter  (subscript  IT  to  assume.  We  then  allow  one  of  the  parameters  to  vary  and 
calculate  the  actual  and  optimal  covariance  with  the  equations  of  Section  !l.  We  define  a 
future  of  merit 


k»i.r 

IV* 

0(t„ . t » = v > 


L +tP/*k> 


(4.17) 


a\  the  integrated  (average)  ratio  of  the  optimal  to  the  actual  predicted  position  error  standard 
deviation.  Since  ax  \<  j is  always  greater  than  or  equal  to  o,  opl.  0 is  always  less  than  unity 
and  equal  to  one  only  for  the  optimal  case.  We  might  therefore  think  of  0 as  the  degree  to 
which  the  siiboptimal  filter  matches  the  optimal  perforr.  •nee.  Using  prediction  time  tp  = 

10  seconds  and  total  run  time  T = 50  seconds  as  usual,  the  sensitivity  results  are  plotted  in 
Figure  4.5  for  three  different  filter  bandwidths  wide  (w).  medium  (M)  and  narrow  (N).  In 
Figure  4.5(a).  we  find  it  takes  an  error  of  approximately  an  order  of  magnitude  in  assumed 
target  maneuver  level  to  produce  a 50  percent  optimality  level  i.e.  when  the  actual  error  is 
twice  what  it  could  be.  It  is  also  slightly  preferable  if  we  must  be  in  error  -to  underestimate 
than  to  overestimate  the  maneuver  level  Oy . In  Figure  4.5(h).  we  find  that  prediction  error 
levels  are  practically  insensitive  to  tire  choice  of  the  maneuver  time  constant  rM . The  reason 
for  (hts  behavior  is  probably  related  to  the  filter-predictor  offsetting  effects  mentioned  in 
the  last  paragraph.  Notice  in  both  figures,  also,  that  there  is  very  little  (if  any)  discernible 
difference  in  sensitivity  our  the  significant  range  of  filter  bandwidths. 

The  author  has  therefore  chosen  the  random-acceleration  target  model  over  the  poly- 
nomial mode!  primarily  Ireeause  the  random-acceleration  model  directly  relates  the  target’s 
kinematics  to  the  filter  transition  and  process  noise  matrices  and  therefore  to  the  filter 
bandwidth.  AH  the  extrapolation  equations  will  change  from  those  presented  for  the  poly- 
nomial case.  Tire  update  equations  do  not  change.  Specifying  the  one-pass  initialization 
completes  the  description  of  the  new  filter.  For  reference  purposes,  a summary  of  tfiesc 
filter  equations  appears  in  Table  4. 1 . 


Before  continuing  to  the  next  section,  it  should  be  mentioned  that  other  nondeter- 
ministic  target  models  are  under  consideration  for  the  GFCS  tracking  application.  Brown 
and  Price  ( I‘>74).  for  example,  tried  a higher-order  analogy  of  the  model  in  this  section. 
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Figure  4.3.  Random  Acceleration  Figure  4.4.  Random  Acceleration  Filter 

Filter  Steady  State  with  Various  Values  of  om 

Values  as  Function  and  tm  Normalized  Standard 

of  oyj  Deviation  at  Predicted  (10 

second)  Position  Error 


i c . acceleration  rare  is  a exponentially  correlated  random  variable,  but  found  that  it  did  not 
work  as  well  as  the  model  discussed  here.  Moose  (1972)  discusses  a very  interesting  model, 
known  as  a semi-Markov  process,  which  he  applied  to  the  case  of  a maneuvering  submarine. 
The  author  would  like  to  investigate  the  applicability  of  this  model  to  the  air-target  track- 
ing problem. 


Table  4. 1 Random  Acceleration  Filter  Equation^.  A Summary 

(1) 

(2) 

(3) 
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EXTRAPOLATION 

Xj  ik/k  - 1)  = ,\|(k  - 1/k  - 1)  + Atx2  (k.  - 1/k  - 1 ) -r  axj(k  — l/k  — 1 ) 
x2(k/k  - 1)  - x2(k  — 1/k  - 1 ) + pk3  (k  - 1/k  — 1) 
x3(k/k  - 1)  ~ 7x3  (k  - 1/k  - 1) 


I 


Table  4. 1 . Random  Acceptation  Filter  Equations:  A Summary-(Continueil) 


Pn(k/k  - 1) 

P,2(k/k-  1) 

P, 3 (k/k  - I) 

P22Ck/k  - I) 

P23(k/k  - 1 ) 
P33(k/k  - I ) 

GAINS 
Kj  (k)  = P 1 1 
K2(k)  = Pj; 
K3(k)  = P|; 

RESIDUAL 
i>(k/k~  1)  = 


= P,  | (k  — l/k  - 1)  + 2AtP,,(k  - !/k  - I ) + 2aPl3(k-  1/k-  1) 

+ At2P22(k-  l/k  - 1)  + 2crAtP23(k-  l/k  - 1 ) 

+ a2P33(k  - l/k  - 1)  + Q,,(k-  1) 

= Pi2(k  - l/k  - 1)  + /JP13(k-  l/k  - 1)  + AtP22(k-  l/k  - 1) 

+ (a  f0At»P,,(k-  l/k  - 1)  + 20P33(k-  l/k  - 1)  + Ql2(k  - 1) 

= -yP,3(k  - l/k  - 1)  + 7AtP23(k-  l/k  - I)  + £*7P33(k  - l/k  - I) 
+ Q,  3(k  - 1) 

= P22(k - l/k  - 1)  + 20P23(k-  l/k  - I)  + 02P33(k-  l/k  - I) 

+ 022^  — ^ ) 

= 7P23(k-  l/k"  I)  + 07P33(k-  l/k  - 1)  + Q23(k  - 1) 

- 72P33(k  - l/k  - I)  + Q33(k  - I) 

(k/k  - l)/[P,,(k/k  - 1)  + o2(k)| 

(k/k  - l)/(Pn(k/k-  1)  + a 2 ( k ) 1 
(k/k  - D/IP, ,(k/k - I)  + o2(k)l 


z(k)  - X|  (k/k  - 1) 
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(4) 

(5) 

(6) 

(7) 

(8) 
(9) 

(10) 

(ID 

(12) 

(13) 


1 
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Table  4.1.  Random  Acceleration  Filter  Equations:  A Summary~(Continued) 
UPDATE 

x,(k/k)  = x,  (k/k  - 1)  + K,(kMk/k-l) 
x2 (k/k)  = x2 (k/k - 1)  <-  K2(kMk/k-l) 
x3(k/k)  = x3 (k/k  - 1)  K3(k)J>(k/k-  1) 

P,  j (k/k)  = (1  — K(  <k)l  Pj  i (k/k  — 1) 

P12(k/k)  = [ 1 - Kj  (k)) P, 2(k/k  - 1) 

PI3(k/k)  = ( 1 - K,(k))Pj3(k/k-  1) 

P22(k/k)  = P22(k/k-  1)  - K2(k)P,2(k/k-  1) 

P23(k/k)  = P23(k/k-  1)  - K2(k)P13(k/k-  1) 


(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

(20) 
(21) 


P33(k/k)  = P33(k/k  - 1)  - K3(k)P,3(k/k-  1) 


INITIALIZATION 
X|  (0/0)  = z(0) 

x2(0/0)  = x3(0/0)  = 0 

P,,(0/0)  = o2(  0) 

P,,(0/0)  =-■  0. 

P13(0/0)  = 0. 

n22(0/0)  = 01(0) 


(23) 

(24) 

(25) 

(26) 

(27) 

(28) 
(29) 


P23(0/0)  = 0. 
P33(0/0)  = 0^(0) 


(30) 


CONSTANTS 
1 = exp  (~At/rM  ) 

0 = rM<  ! -7) 
a ~ T5i(7+  At :/rM  - I) 

PROCESS  NOISE  (CONSTANT) 

= -AtOy /Ty 

Q23  = AtQ„/2 
022  = 2At02J/3 
0,3  = Q;:/2 
0,2  = 3AtQj  j/4 
0„  = 2AI0.2/5 

PREDICTION 

Xj(t  + tp/t)  = X(  (t/t)  + X2<t/t)‘t,,  + X3  (t/t)Ty  |exp  (-tp/Ty  )' 

+ *p/rM  " !l 

RECOMMENDED  VALUES  OF  PARAMETERS 


16  Hertz 

(41)  | 

i 

= 0. 1-10.0  yards/second2 

(42)  j 

f 

= 3-20  seconds 

(43)  | 
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Table  4. 1 . Random  Acceleration  Filter  Equations:  A Summary-(Continued) 


FOR  AIR  TARGETS 


o2(0)  = 300  yards/second 
03 (0)  = 10  vards/second2 


FOR  SURFACE  TARGETS 


o2(0)  = 14  yards/second 
o3(0)  = 2 yard/second2 


V.  ADAPTATION 


The  Kalman  filter  formulation,  presented  in  the  previous  sections,  assumes  complete 
knowledge  of  the  linear  dynamic  model  and  the  process  noise  covariance.  In  a general 
tracking  filter  application,  such  as  the  gunfire  control  problem,  the  particular  strategy  being 
exercised  by  the  target  is  unknown.  The  form  of  the  state  vector  and  its  propagation 
characteristics  is  assumed  and  may  or  may  not  adequately  represent  the  true  target  motion 
over  long  periods  of  time.  (We  always  expect,  however,  the  dynamics  model  to  be  a good 
approximation  of  target  motion  over  short  periods  of  time.)  Such  a situation  is  usually 
referred  to  a ”suboptima!  modeling”  in  the  sense  that  no  attempt  is  made  to  fully  model  the 
target  dynamics.  The  utilization  of  a suboptimal  model  often  leads  to  large  estimation 
errors -a  condition  known  as  filter  divergence.  When  divergence  occurs,  an  inconsistency 
between  the  error  covariance  calculated  by  the  filter  and  the  actual  error  covariance  occurs. 
Examples  of  such  divergence  problems  will  be  shown  shortly. 


An  adaptive  filter  is  basically  a method  of  adjusting  parameters  in  order  to  effect  a 
more  realistic  match  between  the  calculated  and  actual  filter  error  covariances.  The  purpose 
of  such  a technique  is  to  reduce  and  bound  the  actual  error  covariance  when  modeling  errors 
become  large  enough  to  seriously  affect  the  performance.  We  wilt  find  that,  contrary  to  the 
linear  Kalman  filter  techniques  described  thus  far.  the  calculated  error  covariance  must  be- 
come a function  of  the  actual  data  through  a coupling  of  the  filter  parameters  with  the 
target  motion.  Wc  will  also  find  that  the  performance  of  a particular  adaptive  filter  is  a 
function  of  the  procedure  used  to  detect  divergence  and  of  the  method  used  to  modify  the 
filter  when  divergence  is  detected.  After  the  following  examples,  techniques  of  detecting 
divergence  and  two  forms  of  adaptive  filters  will  be  presented. 
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A.  DIVERGENCE  AND  THE  BANDWIDTH  TRADEOFF 


We  will  now  consider  a few  particular  examples  of  target  motion-actually  driven  by 
target  acceleration  that  can  lead  either  to  a filter  divergence  problem  or  perhaps  to  an  un- 
acceptably wide-banded  filter  implementation.  Using  these  examples,  we  will  then  demon- 
strate and  discuss  several  possible  techniques  of  filter  adaptation  to  deal  with  the  situation. 

Let  us  first  consider  the  random-acceleration  filter  model  developed  in  Section  IV.  In 
that  section  we  discussed  the  values  of  the  parameters  representing  target  maneuverability 
and  maneuver  frequency  and  found,  as  might  be  expected,  that  no  single  set  of  numbers 
would  adequately  represent  the  target  scenario.  Rather,  it  was  found  that  a range  of  each 
parameter  could  be  expected  and  that  we  can  essentially  bound  the  parameters  by  a low 
maneuver  level,  low-frequency  (long  time  constant)  parameter  set  which  we  will  designate 
(oma  • w\i,\  ) (where  is  l/Vy  ) and  a high-maneuver-level,  high-frequency  set  (omb  . «mb  >• 
Our  fundamental  assumption  in  approaching  adaptation  with  the  random-acceleration  model 
is  therefore 


(o.MA  • ^ (omb  • <o.mb  1 


where  <om  . cjm  ) is  the  parameter  set  representing  any  actual  target.  Recall  that  a non- 
zero choice  of  (Oma  . coma  * is  due  to  the  fact  that  an  actual  target  moving  through  an  actual 
atmosphere  will  be  buffeted  by  turbulence  (and  perhaps  other  effects)  so  that  (<jm,v  . «*>ma  > 
really  might  represent  inadvertent  maneuvers  of  the  target.  It  is  important  that  we  do  not 
allow  the  A parameters  to  vanish  since  steady-state  error  covariance  would  also  vanish  and 
the  filter  could  diverge  due  even  to  the  mildest  maneuvers. 

Let  us  now  construct  filters  based  upon  each  hounding  parameter  set  and  called  the 
corresponding  (fixed-parameter)  Kalman  filters  A and  B.  We  will  now  exercise  these  filters 
against  several  target  profiles  to  assess  their  acceleration  tracking  performance.  In  Figure  5.1 , 
the  actual  acceleration  of  several  targets  is  shown  by  the  solid  lines.  The  asterisks  plot  the 
acceleration  estimates  of  the  A filter  and  the  circles  those  of  the  B filter.  The  maneuver 
parameter  sets  used  for  these  filters  arc  (o.ma  • wma  I of  (0.5.  I '20)  and  («m  b . u:\t  b > of 
(5.0.  M0».  These  acceleration  profiles  and  their  two  integrals  were  generated  on  the  com- 
puter. and  numerically  generated  white  noise  \.dh  o ~ 5 yards  was  added  to  the  true  posi- 
tions to  simulate  the  measurements.  It  is  observed  that  the  A filter  is  narrow-handed  in  that 
the  acceleration  estimates  are  relatively  smooth.  Unfortunately,  the  A filler  values  tend  to 
diverge  (or  al  least  lag  severely)  from  the  true  value  whenever  the  acceleration  changes 
rapidly.  The  B filter  estimates  tend  to  usually  be  unbiased  (since  filter  B is  very  wide-banded) 
and  never  really  divergent.  Unfortunately,  the  B filter  estimates  always  contain  a large 
amount  of  noise,  even  when  we  see  that  acceleration  can  be  tracked  smoothly  and  accurately 
by  the  A filter.  Figure  5.1  displays  dramatically  the  classical  problem  of  determining  the 
proper  filter  bandwidth  to  yield  the  smoothest  unbiased  estimates.  The  paradox  is  that. 
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while  each  filter  has  its  advantages  and  disadvantages,  neither  is  really  suitable  (as  they  are) 
for  implementation  as  a GFCS  tracking  filter,  it  is  therefore  clear  that  we  should  look  for 
some  method  of  adapting  on-line  the  bandwidth  of  the  filter  to  effect  the  desired  performance. 

Before  leaving  Figure  5.1.  the  reader  might  note  that  these  particular  acceleration 
profiles  obviously  belie  our  target  model  assumption  that  acceleration  is  a stationary,  first' 
order  random  process.  Indeed,  many  real  targets  will  not  follow  this  model  any  more  than 
they  might  follow  a polynomial  model.  The  point  again  is  that  we  want  each  trajectory  to 
be  reasonably  represented  by  the  statistics  of  the  random  process. 


B.  RESIDUAL  STATISTICS  AND  MANEUVER  DETECTION 

When  divergence  occurs,  the  error  vector  (estimated  state  minus  true  state)  grows  large. 
Fortunately,  we  are  able  to  monitor  at  least  partially  the  actual  performance  of  the  filter 
at  any  given  time.  Tins  is  done  by  observing  the  sequence  of  residuals  (often  referred  to  as 
the  “innovations  sequence”)  and  attempting  to  detect  the  buildup  of  a bias  and  consequent 
growth  of  the  residuals.  We  can  determine  the  statistics  of  the  residuals  by  recalling  the 
definition 


£tk/k-  I)  = /(k)  - Hik)  x(k/k-  I) 
and  the  measurement  model 


/ik)  = H(k>  x<k)  ♦ v<k» 

Substituting  for  /( k ).  we  find 

#k/k-  H = wk>  - IKk) £<fc/k  - !)  (5.2) 

By  taking  expected  values,  we  rmd  immediately  that 

I:  fwk/k  - h|  = U (5J) 


and 


i;  iwK/k  hi>f(k/k-  hi  = :«k) 


< 5.4) 


+ Ilfk)  Pik/k  - DlD(k) 

Substituting  the  matrices  for  our  model,  we  find  that,  when  the  (liter  is  operating  optimally, 
the  residual  sequence  should  be  zero-mean  Gaussian  with  variance 
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Op  (k)  = oH k)  + Pi , (k/k  - I) 


(5.5) 


It  is  then  a relatively  simple  matter  to  determine  the  probability  that  the  sampled  residual 
belongs  to  the  population  with  the  above  statistics. 

An  important  extension  of  Equation  (5.2)  can  be  made  in  the  case  of  the  actual 
sampled  residual.  The  actual  residual  is  comprised  of  exactly  the  same  error  terms,  am*  the 
sample  expected  value  (denoted  E-)  is  related  in  the  same  manner  to  the  actU’d  error  <nd 
measurement  covariance,  t.c.. 

f:\  Irtk/k  - \)v}( k/k  - l)|  = Rac,  (k)  + H4rr(k)  PArrfk/k  ~ DHjCT(k)  (5.6) 
or  for  our  ease 

E,  l^tk'k  - l)j  = o;ACr(k> 

(5.7) 

- 'r*A(  T(k)  + pi»A<  r<k/k  ■ S) 

If  we  are  sufficiently  confident  in  our  estimate  of  the  measurement  error  variance  to  assume 
that  our  filter  estimate  n-ik‘  equals  oA<-r(fc».  their  we  find  that  by  using  Equations  (5.5) 
and  (5.7)  we  might  make  certain  important  inferences  about  the  validity  of  our  calculated 
error  covariance  and.  if  necessary,  adjust  it  accordingly.  If  our  estimate  of  o\<-T  (k)  is  very 
inaccurate,  then  adaptive  techniques  based  upon  the  use  of  such  information  would  be  in- 
advisable. Another  technique  to  avoid  this  problem  will  be  discussed  laicr  but  was  not 
implemented. 

Rather  than  work  with  ;;n  individual  residual,  greater  statistical  significance  can  be 
obtained  by  considering  several  data  samples.  We  will  use  a sample  mean  which,  defined  for 
an>  variable  f.  is 


1 

k > = m £ Uk) 


It  is  much  more  convenient  to  implement  the  sample  mean  recursively  by  means  of  a fixed 
memory  length  averager  that  will  approximate  the  exact  sample  mean  to  QtAU.  This  re- 
cursive sample  mean  is  given  as 

ftk)  = Gj  t'tk - I)  + G2  flk)  t>9 


where  the  constant  gams  arc- 
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We  find  that  M • is  the  chi-squared  variable  with  M degrees  of  freedom.  The  expected 
value  of is  unity  and  the  variance  is  two.  It  is  therefore  a relatively  easy  problem  to  con- 
struct tests  concerning  these  variables.  Other  maneuver  detectors  can  be  constructed  by 
consiv'ering  the  correlation  of  residuals.  For  example,  it  can  be  shown  that  the  autocovariance 

flp(i)  = K |*nk/k-  l»r,T(k-  i/k-  I - ill 

sltould  vanish  for  i =£  0.  This  information  forms  the  basis  for  a slightly  different  type  of 
maneuver  detection.  For  example,  the  adaptive  filter  in  the  MARK  86  GFCS  monitors  the 
signs  of  the  residuals  and  declares  a maneuver  when  a certain  number  of  successive  residuals 
show  the  same  sign.  Such  methods  work  very  well.  There  are  other  more  complicated  vari- 
ations of  this  type.  AH  of  these  maneuver  detection  statistics  are  similar  in  that  each  starts 
with  the  assumption  that  the  residuals  should  be  uncorrelated,  zero-mean.  Gaussian  when 
the  filter  is  operating  properly.  In  fact,  the  author  could  find  no  particular  advantage  of 
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any  one  maneuver  detector  over  another  as  far  as  performance  is  concerned.  We  have  chosen, 
for  reasons  of  computational  efficiency,  the  normalized  residual  sample  mean  (Equation 
5.1 1 ) to  construct  our  criterion  for  maneuver  detection.  We  can  now  define  the  maneuver 
detector. 


We  will  define  a maneuver  as  any  target  motion  that  causes  the  (liter  performance,  as 
measured  by  . to  exceed  some  specified  value.  Namely,  a maneuver  is  declared  if 


%<k)l  >Cajrs(k)  = CA/CT 


(5.13a) 


or  equivalently 


l£(k)  >C2i M 


(5.13b) 


where  C is  a constant  that  determines  the  significance  of  the  test.  The  probability  of  in- 
dicating a maneuver  when  there  is  none  (a  “false  detection”  or  Type  I error)  is 


Pi  n ={P?*(k>  >C2op*N(k)} 


(5.14) 


Values  of  P|  j>  as  a function  of  C can  be  found  in  most  introductory  statistics  books.  For 
example,  maneuver  detection  at  the  two-  and  three-sigma  levels  (C  = 2 or  3)  yields  Ppp  = 
4.56"  and  0.2/QT  respectively.  In  order  to  choose  a value  of  C.  we  must  consider  the 
false  detection  probability  in  conjunction  with  the  cost -presumably  in  degraded  performance 
of  such  false  detections.  Such  costs  are  3 function  of  »hc  type  of  action  taken  to  adapt 
when  a maneuver  is  declared.  These  costs  will  be  considered,  at  least  qualitatively,  for  the 
various  adaptive  techniques  to  be  discussed.  A tradeoff  is  involved  since  there  is  a cost  in 
increased  maneuver  detection  time  when  the  maneuver  threshold  C is  raised.  It  is  possible 
(o  estimate  analytically  the  functional  dependence  be  twee,  maneuver  detection  time  and 
the  various  parameters  involved  with  detection  for  certain  simple  maneuvers.  For  example, 
let  us  consider  tile  step  acceleration  target  (of  the  type  discussed  previously)  under  the 
assumption  that  the  filter  is  completely  converged  and  perfectly  nonnrsponsive.  i.c..  the 
filter  error  covariance  and  gains  are  zero.  A step  in  acceleration  a,  will  cause  a residual  bias 
buildup  of  magnitude 


^hudtS)  — , 2,  t^ 


(5.15) 


where  t4  is  the  time  since  application  of  a*.  The  bias  of  the  normalized  residual  sample 
mean  will  then  be 
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*JNbi»«t,)=T~  f *“ait2dt  =~^r-  t3  <5  ,6) 

'III  Jg  2 6 o,Tji  * 

Maneuver  detection  occurs  at  time  To  when 

*^bia»fTo)  * C/y^ST  (5.171 

so  that 

Tp  = (6  T#  C o,l y/M  a,)I/3  <5. 18) 

Using  Equation  (5.5)  and  recalling  that  Pj  j is  assumed  essentially  zero,  we  can  write 

o„  = y/o1  + P|  | = a 

Also  so  that 


TM/vAt  = 

To  £ (6 Co 


(5.20i 


Similar  procedures  can  be  used  to  estimate  maneuver  detection  time  for  other  target  models. 
For  example,  an  acceleration  ramp  (jerque  j\>  would  yield  a residual  bias  of 


I’tiujllj)  ~ ^ jl  tj 


(5.21) 


and  we  find 


Tp  = (24C  o At >* »'**  f 5.22 1 

Equations  such  as  (5.20)  and  (5.22)  are  admittedly  not  exact  since  the  steady-state  filter  co- 
variance  and  grins  are  not  zero  and  the  filter’s  response,  however  small,  would  tend  to  in- 
crease the  time  required  for  detection.  This  effect  is  small,  however,  for  a narrow-bandwidth 
filter.  The  results,  on  the  other  hand,  are  quite  interesting.  The  detection  time  for  both 
cases  is  found  to  be  directly  proporationa!  to  (the  fractional  powers  of)  C.  M and  o and 
inversely  proportional  to  the  magnitude  of  the  step  change.  This  is  intuitively  satisfying 
since  wc  would  expect  it  to  take  longer  to  detect  a small  maneuver  than  a large  one  or 
longer  to  detect  a maneuver  with  a higher  threshold,  a longer  memory  length  or  when 
observed  with  the  superposition  of  more  noise.  Since  the  probability  of  a false  detection  is 
a function  only  of  C.  it  is  seen  that  the  residual  average  does  not  improve  our  maneuver 
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detector.  Indeed,  the  residual  mean  serves  only  to  increase  the  maneuver  detection  time 
without  effect  on  the  false  detection  rate.  Therefore,  the  sample  residual  mean  will  be 
selected  with  M = I.  i.e..  no  memory  length.  If  we  had  chosen  another  type  of  maneuver 
detector  (such  as  estimating  a residual  trend),  this  would  not  necessarily  be  true.  Notice 
also  that,  for  the  step  acceleration  target,  an  increase  in  (he  value  of  C.  say  from  I to  3.  in- 
creases the  maneuver  detection  time  less  than  7 percent  but  improves  the  false  detection 
rate  by  almost  a factor  of  1 7. 

We  have  chosen  liquation  t5.i3)  as  our  maneuver  detector.  It  can  be  shown  to  yield 
t performance  equivalent  to  any  of  the  other  methods  discussed  in  this  section.  Results  with 

this  detector  will  be  presented  in  the  following  sections. 


C.  VARIABLE  BANDWIDTH  ADAPTATION 

Once  divergence  (or  a "maneuver."  as  we  call  any  target  motion  that  produces  divergence 
even  temporarily)  has  been  detected,  a method  of  modifying  or  adapting  the  filter  parameters 
to  correct  the  situation  must  be  specified.  There  are  several  methods  of  dealing  with 
divergence,  most  of  which  effectively  increase  the  gains  to  make  the  filter  more  sensitive  to 
new  data  ar.d.  of  course,  more  sensitive  to  noise.  W'e  will  consider  a brief  survey  of  these 
various  techniques  (excluding  parallel  filtering  to  be  discussed  in  the  next  section). 

Several  survey  and  comparison  papers  on  adaptation  appear  in  the  literature.  Mehra 
( 1972)  classified  the  different  methods  into  four  categories  and  discussed  the  relationships 
between  them  and  the  difficulties  associated  with  each.  Ilagar  1 1973)  conducted  a fairly 
comprehensive  investigation  of  the  various  types  of  adaptive  algorithms  and  their  capabilities 
in  a very  useful  reference.  Sidar  anti  Bar-Shlomo  1 1972)  simulated  and  compared  a number 
of  adaptive  filters  for  appl?  -ation  to  their  gyro  compass  problem.  Particular  emphasis  was 
placed  on  the  Ja/.winski-tyj  • of  covariance  matching  (to  be  discussed  shortly). 

There  arc  several  variants  of  a technique  which  the  author  refers  to  as  bias  correctors. 
Demetrv  and  Titus  1 1 9(»8)  suggested  a bias  corrector,  whereby  if  a maneuver  is  detected  one 
reprocesses  the  most  recent  data  with  a widcr-bandwidlh  filter.  Fricdland  ( 1909)  devised  a 
technique,  whereby  the  state  is  augmented  with  a bias  vector  which  is  then  estimated  in  an 
effectively  decoupled  estimator.  McAulay  and  Denlinecr  ( 1973)  suggested  a multiple-order 
derivative  polynomial  filter  whereby  the  lower  order  (n  = I ) would  be  used  unless  a maneuver 
is  detected.  If  a maneuver  is  detected,  then  a higher-order  (n  = 2)  filter  is  initialized.  Es- 
sentially ihis  technique  attempts  to  break  (he  trajectory  into  piecewise  polynomials  with 
the  break-points  defined  on  line. 

Other  investigators  have  tried  to  directly  adapt  the  memory  length  or  the  gain  matrix 
itself.  Epstein  C 1 97 1 > used  the  memory  length  of  the  filter  as  the  adaptive  variable  and 
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varied  the  memory  as  a function  of  the  residual  series.  Mehra  (1972)  presented  a new 
algorithm  for  the  direct  estimation  of  the  optimal  gain.  Hampton  and  Cooke  (1973)  also 
designed  an  adaptive  filter  for  tracking  high-performance  maneuvering  targets.  This  tech- 
nique uses  the  orthogonality  property  of  the  residual  sequence  to  automatically  track  the 
optimal  gain  levels. 

There  are  some  interesting  algorithms  that  attempt  to  “learn”  the  dynamic  they 
process  the  data.  Mehra  ( 1 97 1 , # 1 ) devised  a technique  to  actually  estimate  on  lute  es- 
sentially all  of  the  (assumed)  linear  system,  i.e.,  the  order,  the  transition  matrix,  and  the 
measurement  and  process-noise  matrices,  even  when  the  processes  are  nonstationary  (of  a 
certain  type).  The  author  has  not  actually  tried  this  algorithm  but  strongly  suspects  there 
would  not  be  enough  time  or  computational  capacity  to  use  this  technique.  Of  course,  one 
could  i.ever  be  certain  that,  once  a strategy  based  on  past  data  was  determined,  the  strategy 
would  continue  to  be  employed  in  the  future.  But,  of  course,  that  factor  is  merely  a part 
of  the  difficulty  faced  by  the  predictor  designer  in  any  case. 

Divergence  can  be  considered  as  having  been  caused  by  an  inaccurate  estimate  of  the 
noise  covariances.  It  is  only  natural  then  that  attempts  are  made  to  estimate  these  covariances 
as  well  as  the  state.  Mehra  (1970)  introduced  a method  to  simultaneously  estimate  Q and  R 
when  the  state  model  was  assumed.  V/eiss  ( 1 970)  surveyed  and  discussed  techniques  of  this 
same  type.  The  principal  objective  of  these  algorithms  is  to  effect  a correspondence  between 
the  actual  covariance  and  the  calculated  covariance-  hence  the  name  for  these  methods, 
“covariance  matching."  Nahi  (1972)  and  Soeda  and  Yoshimura  (1973)  developed  pro- 
cedures to  more  or  less  -optimally  modify  the  calculated  error  covariance  to  prevent  diver- 
gence when  the  residual  is  not  likely  to  have  come  from  the  calculated  distribution.  A 
pioneer  of  the  covariance  matching  technique,  Jazwinski  ( 1 969),  developed  the  concept 
(usually  known  by  his  name)  of  preventing  divergence  by  covering  modeling  errors  with 
noise  and  adaptively  estimating  the  noise  level.  The  Jazwinski  method,  at  least  in  a con- 
cept* *1  level,  is  the  approach  selected  for  our  application.  It  is  particularly  well  suited  for 
our  purpose  since  the  random  acceleration  model  tells  us  much  about  the  proper  structure 
of  the  noise  process  and  its  relation  to  the  target  kinematics. 

The  Jazwinski  technique  uses  “re«*dual  feedback"  to  specify  the  proper  level  of  process 
noise  to  add.  The  process  works  as  follows.  In  the  event  of  a maneuver  detection,  wc  make 
use  of  Equation  (5.6)  (assuming  accurate  knowledge  of  the  measurement  noise  covariance) 
to  estimate  the  level  of  actual  process  noise.  One  assumes  that  the  actual  extrapolated  error 
covariance  consists  of  three  ♦erms,  i.e., 

P,UT(k/k-  1)  = 0(k,  k - l)P(k-  I/k-  1) 0T(k  k-  I)  + Q<k-  1)  + Q*(k- 1 ) (5.23) 

where  the  first  two  terms  on  the  righthand  side  represent  the  calculated  error  covariance 
with  a small  process  noise  covariance  Q assumed  by  the  filter.  The  term  Q’  is  effectively 
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added  to  balance  the  equation.  If  there  is  no  maneuver  detected,  Q*  is  zero.  It  should  be 
noted  that  it  is  not  possible  to  uniquely  specify  Q*  with  this  equation.  We  can  only  deter- 
mine HQ*  HT  orQjj  for  our  case.  Tliis  is  not  particularly  troublesome  as  we  can  choose 
some  normalized  form  of  the  0*  matrix  such  as  the  random  acceleration  model  process 
noise  matrix  to  automatically  define  the  other  elements  of  Q*  as  a function  of  Qj , . For 
our  case,  we  could  calculate 

Qu(k-  1)  = (7$(k)-  1/Ml  (o2(k)  + Ptt(k/k-  1)1  (5.24) 

when  a maneuver  « s detected,  calculate  the  remaining  elements  of  0*  with  the  Singer 
form,  and  add  Q’ » he  error  covariance.  This  technique  was  simulated  and  was  found  to 

yield  fairly  good  -five  filter  performance.  Unfortunately,  this  adaptive  filter  would  oc- 

casionally display  erratic  behavior  and  generate  unreasonable  error  covariances.  In  order  to 
eliminate  this  problem  and  properly  constrain  the  error  covariance,  a different  means  of 
specifying  Q*  was  selected. 

We  therefore  return  to  our  original  assumption  that  the  target  maneuver  leve!--and  sub- 
sequent process  noise -is  bounded  by  the  A and  B parameter  sets.  Therefore,  in  the  absence 
of  maneuver  detection,  process  noise  corresponding  to  set  A will  be  added  and.  if  a maneuver 
is  declared.  Qlt . corresponding  to  set  B.  will  be  added. 

That  is: 


if  [?*  <k)  < C2/M I . then  Q(k  - 1 ) = Qa  . (5.25a) 

if|^(k)  > C-/M),  then  Q(k  - 1)  = Ob  (5.25b) 

If  a maneuver  is  declared,  some  time  will  elapse  before  QB  builds  up  the  error  covariance 
matrix.  Another  obvious  alternative  is  to  simply  reset  the  error  covariance  in  order  to  obtain 
a faster  l espouse.  This  technique  was  rejected,  however,  since  one  must  pay  a very  large 
price  for  a false  detection.  Using  Liquation  (5.25).  it  is  entirely  possible  for  the  maneuver 
WU\  or  to  turn  off  before  the  error  covariance  builds  up  to  the  steady-state  value.  Con- 
e i .ally.  Figure  5.2  shows  the  sequence  of  events.  Initially,  the  target  is  not  maneuvering 
unu  the  filter,  using  the  A parameters,  is  tracking  acceleration  very  well.  After  the  maneuver 
occurs,  an  interval  of  time  elapses  before  the  maneuver  is  detected  which  we  have  been  re- 
ferring to  as  Tdtft«cth>R  - Once  the  maneuver  is  detected  and  we  start  adding  Ob  • >t  takes 
another  amount  of  time,  called  Tadapt.  for  the  filter  to  "open”  or  to  increase  the  error  co- 
variance  sufficiently  to  remove  the  bias.  When  the  maneuver  detector  decides  the  bias  is 
gone  and  turns  off.  the  filter  returns  to  adding  the  small  process  noise  Qa  . The  time  re- 
quired to  return  to  steady  state  with  the  A parameters  is  called  Tfcgonm(:c. 
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Figure  5.2.  Conceptual  Example:  Single  Variable-Band  Adaption 


In  Figure  5.3.  the  example  acceleration  profiles  are  shown  being  tracked  by  the  single, 
adaptive,  variable-bandwidth  filter.  The  maneuver  parameter  sets  are  the  same  as  those  used 
in  the  previous  nonadaptive  examples.  The  maneuver  detector  was  operated  at  the  ('  = 3 
(sigma)  level  to  minimize  false  detections  and  no  memory  (M  = I ) to  minimize  detection 
time.  It  is  observed  that  this  technique  tracks  rather  smoothly  while  improving  the  bias 
error.  There  tends  to  be  some  overshoot  when  the  filter  adapts  because  it  was  initially 
lagging  and  builds  up  excessive  rates  in  order  to  catch  up.  The  filter  appears  to  respond  well 
(as  expected)  to  the  step  change  in  acceleration.  The  predicted  detection  time  of  1.65  sec. 
apparently  matches  the  simulation  quite  well.  The  variable-bandwidth  or  adaptive  filter 
represents  a marked  improvement  in  accuracy  over  either  of  the  fixed-bandwidth  filters 
discussed  previously. 


D.  DUAL  BANDWIDTH  ADAPTATION 


Upon  considering  the  single-variable-bandwidth  filter  just  discussed,  it  was  felt  that 
certain  improvements  could  be  made  in  the  adaptive  performance.  Ideally,  when  a maneuver 
has  been  detected,  one  would  likely  reprocess  the  measurements  over  some  interval 
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Figure  5.3.  Acceleration  Tracking  Performance:  Examples  with  Adaptive 
Variable  Bandwidth  Filter 
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immediately  preceding  current  time  with  a wider-bandwidth  filter  so  as  to  remove  the  bias 
error  which  presumably  has  been  occurring.  Unfortunately,  our  application  essentially  de- 
mands a fully  recursive  filter  in  order  to  efficiently  implement  the  same  on  a limited  real- 
time computer.  It  would  be  very  difficult  to  interrupt  normal  processing  in  order  to  re- 
process past  data.  The  best  answer  obviously  is  some  type  of  recursive  wide-bandwidth 
filter  operating  in  parallel  to  the  “main  filter"  which  can  be  used  to  (more  or  less)  instan- 
taneously remove  the  bias  error  once  it  has  been  detected.  Such  a technique  has  the  obvious 
advantage  of  eliminating  the  last  two  waiting  intervals  of  the  single  variable  bandwidth 
filter- the  adaptation  and  reconvergence  times.  If  done  properly,  it  can  also  eliminate  the 
overshoot  after  adaptation. 

The  dual-bandwidth  filter  would  work  as  follows.  Two  filters.  A and  B.  corresponding 
to  the  respective  maneuver  parameter  bounds,  would  operate  simultaneously.  The  filter 
would  only  output  (to  the  FCS)  the  state  vector  of  filter  A.  Xa  • If  divergence  of  filter  A 
is  detected,  using  the  detection  criterion  of  Equation  (5.25).  the  state  vector  of  filter  B. 
xb  • which  should  be  unbiased,  is  put  into  filter  A.  i.e.. 


Iflt^Afk)  > C2 /M 1 . then  xa  - xb 


(5.26) 


Conceptually,  we  want  the  adaptation  to  work  as  in  Figure  5.4.  The  A filter  is  outputting 
smooth  estimates  of  the  state  until  the  maneuver  is  detected.  Using  Equation  (5.26).  the  A 
filter  then  “jumps”  to  the  current  (unbiased)  estimate  of  the  B filter  and  no  adaption  or 
reconvergence  is  required.  Admittedly,  the  output  vector  xa  will  be  discontinuous,  but  in 
this  situation  xb  represents  the  best  information  available.  An  important  consideration 
was  the  decision  as  to  what  modification,  if  any.  should  be  made  to  the  A filter  bandwidth. 
Three  options  were  considered  and  tested. 

Theoretically,  if  one  resets  the  state  estimate  of  A to  that  of  B.  the  bandwidth  should 
be  similarly  reset.  That  is.  if  a maneuver  is  declared  and  Equation  (5.26)  is  in  effect,  then 


Option  I:  PA  = P« 

This  option  was  found  unacceptable,  however,  as  one  pays  a high  cost  of  a false  detection 
since  the  long  reconvergence  time  has  not  been  eliminated.  Leaving  the  A bandwidth 
unchanged,  i.e.. 


Option  2:  PA  " Pa 


was  desirable  since  it  did  net  suffer  the  disadvantages  of  Option  1 . The  particular  maneuver 
detector  we  are  using,  however,  did  not  do  well  with  Option  2 since  the  large  random  errors 
of  the  B filter  look  like  biases  to  the  A filter  and  repeated  maneuver  detections  tended  to 
occur.  In  order  to  eliminate  this  problem,  the  author  decided  on 
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Figure  5.4.  Conceptual  Example:  Dual-Bandwidth  Adaptation 

Option  3:  PA  = PA  + Ob  (5.27) 

i.e..  a gradual  widening  of  the  A Filter  bandwidth  which  eliminated  the  cyclic  detections  and 
maintained  a low  cost  of  a false  detection.  In  Figure  5.5,  the  same  acceleration  profiles 
were  simulated  with  the  same  parameters  for  the  maneuver  statistics  and  maneuver  detector 
but  with  the  dual-bandwidth  filter.  The  acceleration  estimatioi.  errors  were  reduced  even 
further  and  the  overshoot  errors  disappeared.  The  root-mean-square  acceleration  estima- 
tion error  was  reduced  on  the  average  by  25  percent  over  the  single-bandwidth  filter. 

Another  interesting  possibility,  only  superficially  examined  by  the  author  to  date,  in- 
volves the  generation  of  an  output  vector,  x,,.  which  is  always  a linear  combination  of  x* 
and  Xr  : i.e.. 


= w Xa  + ( 1 - W)  Xr 


(5.28) 


where  (he  weighting  factor  W is  a function  of  the  residual  statistics  both  sample  and 
calculated  -of  both  the  A and  B filter.  An  obvious  advantage  of  such  a technique  would  be 
to  construct  a ( more-or-less)  optimal  combination  of  the  A and  B filter  with  continuity  of 
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Figure  5.5.  Acceleration  Tracking  Performance:  Examples  with  Adaptive 
Dual  Bandwidth  Filter 
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the  output  state  vector  simply  a by-product  not  without  merit  of  its  own.  Thorpe  (1973) 
developed  a dual-bandwidth  filter  concept  in  which  target  maneuvers,  assumed  to  be  relatively 
infrequent  events  were  modeled  by  introducing  a binary  random  variable  in  the  target  state 
equation.  He  used  the  likelihood  ratio  for  the  detection  of  a maneuver  to  determine  the 
weight  W.  Brown  and  Price  ( 1974)  studied  the  .ility  of  a bank  of  parallel  filters,  each  with 
a different  bandwidth,  using  the  random  acceleration  model.  A combined  estimate  was 
constructed  on  the  basis  of  a hypothesis  test  of  the  probability  that  each  filter  is  the  correct 
one.  Alspach  (1973)  also  constructs  a bank  of  parallel  filters  and  uses  Bayesian  techniques 
to  estimate  the  optimal  gain.  Moose  (1972)  also  uses  Bayesian  methods  to  determine  the 
relative  weighting  of  the  output  of  a bank  of  Kalman  filters  to  form  the  best  combined 
estimate.  The  author  believes  these  techniques  deserve  further  consideration  for  implementa- 
tion in  a GFCS  application. 

In  Equation  (5.7).  we  related  the  sample  residual  variance  to  the  actua!  measurement 
and  estimation  error  variance  and  mentioned  the  problem  that  if  our  estimate  of  the 
measurement  variance  is  poor  then  maneuver  detection  based  upon  that  equation  would  not 
work  well.  To  illustrate,  let  us  consider  the  sensitivity  of  the  adaptation  process  io  targe 
errors  ir.  cur  estimate  of  measurement  error.  If  we  greatly  underestimate  the  measurement 
error  level,  then  the  sample  mean  residual  will  often  exceed  the  expected  value  due  to 
measurement  error,  and  false  maneuver  detections  will  result.  Unfortunately,  such  false  de- 
tections lead  to  a filter  bandwidth  change  in  just  the  opposite  direction  of  the  proper  adapta- 
tion. Instead  of  increasing  the  gains  and  thereby  weighting  the  measurements  heavier,  we 
should  increase  R.  thereby  decreasing  tic:  gains.  Conversely,  if  we  grossly  overestimate  the 
measurement  error  variance,  then  maneuver  detection  is  delayed  (perhaps  indefinitely)  so 
that  we  do  not  increase  the  gains  to  follow  the  target.  For  the  application  intended  for  this 
filter,  this  problem  should  not  occur,  as  good  estimates  of  the  sensor  statistics  should  be 
available  in  order  to  properly  model  them. 

Another  approach  which  the  author  wishes  to  pursue  is  to  statistically  analyze  directly 
the  residual  sequence  for  each  filter  of  a parallel  filter  bank  without  regard  to  the  calculated 
and/or  assumed  statistics.  In  other  words,  we  simply  look  at  each  residual  sequence  and 
pick  the  smoothest  unbiased  one.  This  technique  could  also  be  incorporated  with  the  weight- 
ing technique  of  Equation  (5.28).  The  author  has  experimented  with  such  a rcsidual- 
analysis/weighttng-factor  approach  but  has  not  yet  determined  a method  which  works  well. 
Brown  and  Price  (1974)  also  tried  a variant  of  this  approach  with  reasonable  success.  The 
principal  advantage  of  such  a techmquc.  if  workable,  would  be  that  it  has  essentially  no 
sensitivity  to  an  incorrect  estimate  of  the  measurement  statistics. 


In  summary,  the  author  has  tried  or  at  least  studied  many  different  adaptation  tech- 
niques but  has  to  date  found  none  that  worxs  any  better  than  the  one  presented.  There 
appear,  however,  to  be  a large  number  of  alternatives  in  the  literature  (some  of  them  men- 
tioned in  this  section)  that  offer  promise  of  improvement  and  which  should  be  investigated 
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further.  The  concept  of  a bank  of  filters  of  various  bandwidths  is  particularly  appealing. 
Sensitivity  to  incorrect  observation  error  statistics  is  probably  the  most  serious  problem  that 
must  be  addressed.  This  problem  will  be  discussed  in  a future  report. 


VI.  SERIALLY  CORRELATED  MEASUREMENT  ERROR 

In  Section  11.  the  conventional  discrete  Kalman  filter  was  described.  An  important  and 
necessary  assumption  for  that  development  was  that  the  measurement  error  be  “white'”  or 
uncorrelated.  In  Section  III.  we  found  that  it  is  desirable  to  process  data  at  the  nighest  rate 
possible  when  the  measurement  errors  are  independent.  Unfortunately,  the  assumption  of 
white  noise  and  the  desirability  of  a high  data  rate  are  often  incompatible  when  using 
measurements  from  real  physical  systems.  For  example,  truly  white  noise  never  actually 
exists  in  an  real  system.  When  sampling  data  from  such  a system,  eventually  one  usually 
finds  the  noise  autocorrelated  (serially  correlated)  when  observed  at  some  sufficiently  high 
data  rate. 

In  this  section,  we  will  address  this  problem.  We  will  present  a common  noise  model 
and  a modification  to  the  conventional  Kalman  filter  that  deals  with  this  model.  W:e  will 
find  that  one  can  quantitatively  assess  the  degradation  in  performance  due  to  autocorrelation 
and  can  perform  simple  sensitivity  ar  jlyses  to  determine  the  effects  of  incorrect  estimates  of 
the  statistical  parameters. 

A.  NOISE  MODEL 

The  measurement  model  for  the  Kalman  filter,  as  presented  in  Section  II.  is 

z(k)  = H(k)  x(k)  + Wk) 

where  it  was  assumed  that 

E Iy(k)l  = 0 
and 

E ly(j)  vT(k)|  = R(k)5jk 

We  will  now  replace  the  assumption  of  the  last  equation  with  the  more  general  assumption 
that  the  measurement  noise  v is  the  output  of  a linear  discrete  dynamic  system  driven  by 
white  noise.  I.e.. 

vfk)  = ♦(k.  k - I ) v(k  - I ) + |(k  - I ) (6.1) 
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where- 


E l|(k)J  = 0 


(6.2) 


and 


E «U)r^Ck)l  = R*(k>6jk  (6.3) 

Sage  and  Melsa  (1971)  point  out,  “Although  this  is  not  the  most  general  form  of  colored 
noise,  it  is  probably  the  most  general  practical  form.  One  often  has  extreme  difficulty  in 
establishing  the  parameters  of  such  a simple  model,  so  that  it  is  hard  to  consider  the  use  of 
any  more  general  form."  It  is  easy  to  determine  that 

R(k)  = Wk.k-  1 ) R(k  - l)*T  (k.k-  1)  + R*(k-  I)  (6.4) 

and  that  the  autocovariance  of  v is 

E ly  (k)yT  (k  - 1)|  = +(k.k-  I ) R(k  - I)  (6.5) 

These  equations  are  useful  in  defining  ♦ and  R* . 

Now  let  us  consider  a particular  type  of  first-order  process  chosen  for  the  measure- 
ment error  model  for  this  application.  An  exponential  autocorrelation  function  is  found  to 
be  both  convenient  and  reasonably  matches  power  spectral  density  information  that  has 
been  estimated  for  candidate  sensors  for  the  MARK  68  GFCS.  The  noise  propagation  equa- 
tion for  this  case  is: 


vtk)  = p(k)  v(k  - 1)  ♦ ow(k)|(k-  1) 


(6.6) 


where 


p(k)  = exp  |-  AtMk}|  (6.7) 

is  the  correlation  coefficient  for  lag  t.  The  standard  deviation  of  the  white  noise  driving 
term  is 


o*(k)  = o(k>v/l  - p:<k) 


(6.8) 


and 


o(k)  = E |v*  (k)j 


(6.9) 
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as  before.  We  will  be  primarily  concerned  with  constant  values  of  r between  0 3nd  0.30 
seconds.  Although  we  will  consider  parametrically  a wider  range  of  values,  note  that  this 
process  is  not.  sirickly  speaking,  assumed  to  be  stationary  or  Gaussian  although,  in  actuality, 
at  is  only  slowly  varying  and  slightly  non-Gaussian. 

A digital  noise  generation  program,  required  for  simulation  purposes,  is  discussed  and 
presented  in  Appendix  A.  The  algorithm,  based  upon  a paper  by  L.  F.  Balas  (1967).  will 
generate  noise  of  the  type  assumed  by  Equations  (6.6)  to  (6.8). 


B.  RESTRUCTURED  KALMAN  FILTER 

In  this  section,  the  details  of  applying  an  algorithm  to  process  data  with  autoconrelated 
measurement  noise  (of  the  type  discussed  in  the  previous  paragraphs)  will  be  presented.  The 
algorithm,  from  Sage  and  Melsa  ( 1971 ).  is  reproduced  in  Table  6.1,  Before  substituting  the 
matrices  for  our  system  into  the  algorithm,  it  was  found  convenient  to  rearrange  the  given 
equations  for  our  purposes.  The  original  form  of  the  algorithm  consists  of  two  processes: 
smoothing  and  filtering.  The  author  found  that  it  is  possible  to  rewrite  these  equations  so 
that  the  algorithm  appears  to  be  of  an  extrapolation-update  form  which,  of  course,  is  the 
method  the  conventional  Kalman  filter  in  Section  H is  written  in.  It  was  found  that  this 
form  requires  less  computation  than  the  original  and  can  be  easily  related  to  the  white  noise 
filter.  It  is  pertinent  to  note  at  this  time  that  the  Sage  and  Melsa  algorithm  was  chosen  over 
the  augmented  state  approach  since  the  former  does  not  require  3n  increase  in  the  dimension 
of  the  state  vector  and  does  not  result  in  ill-conditioned  computations.  The  authors  prin- 
cipally responsible  for  the  algorithm  in  Sage  and  Melsa  are  Bryson  and  Henrikson  ( 1968)  who 
discuss  the  relative  merits  of  the  two  approaches. 

Let  us  now  consider  the  smoother/filtcr  algorithm  in  Table  6.1.  Notice  first  that  the 
residual  for  the  smoothing  and  filtering  equations.  Equations  (8)  and  (9).  are  identical.  We 
will  define  this  new  a priori  residual  as 

n*(k/k-  I)  = z(k)  - *(k.  k-  l)z(k-  I)  - H*(k-  hx(k-  l/k  - I)  (6.10) 
If  we  substitute  the  definition  of  H'tk  - I ) from  Equation  (7),  we  find  we  can  write 

£*(k/k  - ! ) = n ( k/k  - ! ) - *«k.  k - llvtk  - l/k  - I)  (6.1!) 

where 

viklk  - I)  = z(k)  - H(k)xfk/k  - I)  ri>.!2) 

is  the  usual  residual  for  the  white  noise  filter  and 
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Tabic  6. ! . Smoother/Filtcr  Algorithm  for  Autocorrelatsd  Measurement  Noise 


MODEL 

x(k  + I ) = #k  + I . k)  x(k)  + w(k)  ( I ) 

OBSERVATIONS 

z(k)  = H(k)  x(k)  + v(k)  (2) 

v(k  + I)  = +(k+l.k)v(k)  + f(k)  (3) 

STATISTICS 

E(w(k)J  = E|£(k)I  = E[w(j)fT(k)j  = o (4) 

EIw(j>wT<k>I  = Q(kl£jk  (5) 

= R*(k)6jk  (6) 

DEFINITION 

H*(k  - I)  = H(k)p(k,k-  I)  - +<k,k-  I)H(k-  I)  (7) 

SMOOTHER 

x(k  - l/fc)  = x(k  I/k  - I)  + K»(k-  l)(2(k) 


- *fk.k-l)z(k-  I)  - H*(k  — I )x(k  - I/k  ~ l)|  (8) 

FILTER 

x(k/k)  “ $(k,  k - I)  x(k  - I/k)  + Kf(k)|z(k) 

- ♦(k,  k - l)z (k  - I)  - H*(k  - I)x(k  - I/k  - I)|  (9) 
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Tabic  6. 1 . Smoother/'Filtcr  Algorithm  for  Autocorrelated  Measurement  Noise -(Continued) 


SMOOTHER  GAiN 


K«(k  - I)  = P(k  - I/k  - l)H*T(k  - l)[H*(k-  l)Pfk-  1/k  - l)H*T(k-  I) 


+ R*(k  - 1)  + H(k)Q(k  - l)HT(k)| 


FILTER  GAIN 


Kr(k)  = Q(k  - l)HT(k)|H*(k  - !)P(k-  l/k-  l)H*T(k-  !) 


+ R*(k  - I)  + H(k)Q(k  - l)llT(k)| 


T/t-ii  '* 


SMOOTHED  COVARIANCE 


Ptk  - l/k)  = (I  - Kt(k  - 1 )H*(k  - I )|  P(k  - l/k  - ! ) 


II  - K,(k  - I )H*(k - 1)|  + KjK  - I )|R*(k  - 1)  + H(k)Q(k-  l)HT(k)j  Kjfk  ~ h (12) 


FILTERED  COVARIANCE 


Pck/k)  = otk.  k - I)P(k-  I/k)0T(k.  k - I)  + Qtk-  D 

- K.(k)|H*(k  - l)P(k  - l/k  - l)H*T(k-  I)  + R*fk  - l) 

+ H(k)Q(k-  l)HT(k)|KfT(k)  - $(k.  k - DKjk-  l)H(k)Q(k  - I) 

- 0<k  - l)HT(k)Kj(k-  l)^T(k.k  - I) 


Wk-  l/k  - I)  = ilk  - I)  - H(k  - I ) x(k  - l/k  - 1) 


(6.13) 


is  the  a posteriori  residual  from  the  previous  filter  cycle.  We  will  now  define 

»?(k)  = *(k.k  - l)£(k-  l/k  - I) 


which  is  essentially  the  correlated  portion  of  the  old  residual  that  we  want  to  remove. 
Therefore.  Equation  (6. 1 1 ) becomes 


f 


n*(k/k-  I)  = 4k/k  - h - n(k) 


(6.15) 


Since  the  smoother  and  filter  equations  are  of  similar  form,  it  is  easily  seen  that  the  two 
equations  can  be  combined  to  produce  “extrapolation"  and  "update"  equations,  identical 
in  form  to  the  white  noise  filter,  l.e.. 

x(k/k  - i>  = <Mk.  k-  I ) Afk  - l/k-  I)  (6.16) 

and 

x(k/k)  = xtk/k  - I)  + K*(k)£*(k/'k-  1)  (6.17) 

where 

K*(ki  = ptk.k-  I ) K,(k  - I ) + Kf(k)  (6.18) 

is  the  equivalent  gain. 

Now  let  us  consider  the  gain  matrices.  Notice  that  the  inverse  of  the  matrix 
Gtk)  - H*tk-  l)l»(k  l/k-  1 ) H*5 (k  - !)  + R*(k-  I)  + H(k)Qtk-  hI|T(k)  (6.19) 
appears  in  noth  gain  equations,  l.e.. 

K<(k  It  = Pi k I k - I ) H*Ttk  - DC  »(ks  (6.20) 

and 

Kp(k)  f 0(k-  DHTtk)C  Mk)  16.21) 

Therefore 

K*(k)  = [<Xk.k-  i ) Ptk  - l/k  - I ) H*T(k  - i)  + CXk  - !)IIDk)|  G - (k»  (6.22) 
By  defining  the  extrapolated  covariance  as  before,  i.e.. 

Ptk/k-  I)  = ptk.  k-  D Ptk  - |/k-  i >p7  (k.  k - It  + Q(k-  I)  (6.2.D 
wc  can  (hen  write 

K *( k > = IPtk/k-  I ) HT(k)  - A(k»|C  Dk)  (6.24) 
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and 


G(k)  = H(k)  P(k/k  - l)HT(k)  + R*(k-  1)  + B(k)  (6.25) 

where 

A(k)  = 0(k,  k - 1)  P(k  1/k-  I > HT (k - l)'i'T(k,k-  1)  (6.26) 

and 

B(k)  = Wk.k-  1)  H(k  1)  P(k - 1/k  - l)HT(k-  l)^T(k,k-  1) 

- H(k)  0(k,  k - 1)  P(k  - l.k-  l)HT(k-  l)*T(k,k-  1)  (6.27) 

- •‘Kk.k-  1)  H(k  - 1)  P(k  - 1/k  - l)0T(k,k-  l)HT(k) 

Obviously,  if 'P  vanishes,  then  both  A and  B vanish  and  R*  - R.  We  therefore  recover  the 
original  form  of  the  white  noise  gain  equation. 

Now  we  will  consider  the  error  covariance  equations.  The  smoothed  error  covariance 
equation.  Equation  (12).  can  be  rewritten  in  a much  simpler  form  by  expanding  tlr 
term  and  recombining,  making  use  of  the  gain  equation,  i.e., 

P(k-  1/k)  = P(k-  1/k-  I)  + Ks(.,  !)  H*(k  - l)  P'k  - 1/k-  l)H*T(k.  j > kJ  (k  - 1) 

- P(k-  1/k-  l)H*T(k-  l)Kj(k-  1)  Ks(k-  1)  H*(k-  l)P(k-  1/k  - 1) 

+ K,(k-l)|R*(k-  1)  + Htk)Q(k-  l)HT(k)l  (k  - 1) 

= P(k - 1/k-  I)  + K,(k-  l)G(k)Kj(k-  1)  “(k-  1/k-  l)H*T(k-  1)  Kj(k-  1 ) 

- K,(k  ',)H*(k-  1 ) P(k  - 1/k  - 1) 

-=  {!-  Ks(k  - 1 ) H*(k  - 1)1  P(k  - l/k  - !)  (6.28) 

Using  the  smoothing  gain  equation  and  the  i'c  :t  that  P and  G are  symmetric,  we  can  also 
.v.ii*:  this  equation  as 

P(k-  1/k)  = P(k - 1/k  - 1)  - Ks{!; - l)G(k)Ks(k-  1)  (6.29) 

Using  the  filter  gain  equation,  we  nn  also  simplify  the  filtered  covariance  equation. 
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P(k/k)  = otx.k-  I ) P(k  - 1/k)  0T(k,  k - 1)  + Q(k  - 1)  - Q(k  - i)  HT(k)  Kf  (k) 

- 0(k.  k-  1)  Ks(k-  1)  H(k)  CHk  - 1)HT  Q(k-  1)  H^(k)  K[(k-  l)0T(k,k-  1) 

(6.30) 

= 0(k.  k - 1 ) P(k  - 1/k)  0T (k.  k - 1)  + Q(k-  1)  - Q(k  - I)  HT(k)  K*T(k) 

- 0(k.  k - I)  Ks(k  - 1 ) H(k)  Q(k  - l) 

By  substituting  Equation  (6.29)  for  the  smoothed  covariance  and  then  for  the  smoothing 
pin.  we  find 

P(k/k)  = 0(k.  k - 1 ) P(k  - 1/k  - l)0T(k.  k-  1)  + Q(k  - 1) 

- |0{k,  k - 1 ) Ks(k  - 1))  G(k)[0(k.  k-  1)  Ks(k  - 1)]T 

- Q(k-  l)HT(k)  K*T(k)  - 0(K,  k - l)Ks(k-  l)H(k)Q(k-  1)  (6.31) 

= Pt k/k - 1 ) - |K*(k)-Q(k-l)HT(k)G",(k)|  G(k)  |K*(k)-Q(k-  l)HT(k)G-*(k)JT 

Q(k-  !)Hr(k)K*T(k)  - (K*(k)-Q<k-  l)HT(k)G  '(k)l  H(k)Q(k-  1) 

= P(k/k  - 1)  - K*(k)  G(k)  K*T(k) 

which  is  of  a form  identical  to  the  white  noise  filter.  A summary  of  the  equations  for  the 
equivalent  extrapolation/update  algorithm  is  found  in  Table  6.2.  Upon  comparing  this 
algorithm  with  the  white  noise  filter  (Table  2.1 ).  we  find  that  the  correlated  noise  algorithm 
requires  the  additional  calculation  of  A(k).  B(k)and  r?(k).  A FORTRAN  IV  version  of  this 
algorithm,  called  CORKAL.  was  writ.en  by  the  author  and  can  be  found  in  Appendix  A. 

We  now  apply  this  algorithm  to  our  model. 

Our  model  will  now  consist  of  the  current  version  of  the  dynamics  and  prediction,  as 
developed  in  the  previous  sections,  but  we  will  add  the  cor,  luted  noise  model.  The  measure- 
ment noise  transition  matrix  is 


'M.k  - 1)  = 'P  = lp(k)|  (6.32) 

The  standard  measurement  noise  convariancc  matrix  R will  be  replaced  by  the  white  meas- 
urement noise. 


R*(k)  = (oi(k)l 


(6.33) 


60 


Upon  substituting  these  matrices  into  the  algorithm,  we  find  that  the  complicated  matrix 
expressions  reduce  to  rather  simple  algebraic  calculations  for  our  problem.  The  matrix  A is 
(3  X 1 ),  and  B and  ti  are  scalars.  The  algorithm  with  the  indicated  modifications  can  be 
found  again  summarized  in  Table  6.3.  (None  of  the  filter  algorithms  considered  in  this 
section  are  adaptive.  The  adaptive  versions  are  exactly  analogous  to  the  white  noise  models 
in  Section  V.) 


Table  6.2.  Equivalent  Extrapolation/Update  Algorithm 
for  Autocorrelated  Measurement  Noise 

MODEL 

x(k  + 1)  -•  0(k  + l,k)x(k)  + w(k)  (1) 


OBSERVATIONS 

z(k)  = H(k)x(k)  + v(k)  (2) 

v(k+  It  * 4Mk  + l.k)v(k)  + £(k)  (3) 

STATISTICS 

E[w(k)|  = E||(k)|  = E|w(j)£T(k)|  = 0 (4) 

E(w(j)wT(k)|  = 0(k)6jk  (5) 

El|(j)|T(k)|  = R*(k)6jk  (6) 

STATE  EXTRAPOLATION 

x(k/k-  1)  = 0(k,  k - l)x(k-  1/k-  I)  (7) 


COVARIANCE  EXTRAPOLATION 

A(k)  = ^(k,  k - l)P(k  - 1/k  - I)HT(k-  l)*T(k.k-  1) 


(8) 


Table  6.2.  Equivalent  Extrapolation/Update  Algorithm 

for  Autocorrelated  Measurement  Noise-(Continued) 

B(k)  = 'Pfk, k - l>H(k - l)P(k-  1/k-  l)HT(k-  l)*T(k,k-  1) 


- H(k)0(k,  k - l)P(k  - 1/k  - l)HT(k  - l)*T(k,k-l) 

- 'Hk,  k - l)H(k  - 1 )P(k  — 1/k  - l)*T(k,  k - l)HT(k)  (9) 

P(k/k  - 1)  = 0(k.  k - l)P(k  - 1/k  - l)*T(k-  1)  + Q(k  - 1)  (10) 

GAIN 

G(k)  = H(k)P(k/k- l)HT(k)  + R*(k-I)  + B(k)  (11) 

K*(k)  = JP(k/k  - 1 )HT(k)  - A(k)lG_1(k)  (12) 

STATE  UPDATE 

x(k/k)  = x(k/k - 1 ) + K*<k)[z(k)  - H(k)x(k/k-l)  - ^(k)|  (13) 

COVARIANCE  UPDATE 

P(k/k)  = P(k/k  - 1)  - K*(k)G(k)K*T(k)  (14) 

A POSTERIORI  RESIDUAL 

mk)  = 'Hk,  k - l)lz(k)  - H(k  - l)x(k-  Ilk  - I))  (15) 


62 


* 


Table  6.3.  Filter  Equations  for  Correlated  Measurement  Noise:  A Summary 


STATE  EXTRAPOLATION 

x,(k/k-  1)  = x,(k-  1/k  — 1)  + Atx2  (k  - 1/k  - 1)  + ax3(k-  1/k-  1)  (1) 

x2(k/k-l)=  x2(k-  1/k  - 1)  + /3x3  (k  — 1/k  - 1)  (2) 

x3(k/k- 1)  = 7x3  (k- 1/k- 1)  (3) 

CORRELATION  MATRICES 

A,(k)  = p(k)[P,  t(k  - 1 'k  - 1)  + AtP12(k  - 1/k  - I)  + aP,3(k  - 1/k  - 1)]  (4) 

A2(k)  = p(k)[P|  2(k  — 1/k  - 1)  + 0Pl3(k-  1/k  - 1)J  (5) 

A3(k)  = p(k)TP,3(k  ~ 1/k  - 1)  (6) 

B(k)  = p2(k)P,,(k-  1/k  - 1)  - 2A,(k)  (7) 

COVARIANCE  EXTRAPOLATION 

P, , (k/k  - 1 ) = P,,(k-  1/k-  1)  + 2AtP,,(k-  1/k-  I)  + 2aPl3(k  - 1/k  - 1) 

+ At2P22(k-  1/k  - l)  + 2oAtP23(k-  1/k  - 1)  + a2P33(k  - I/k  - I) 

+ Oi  j (k  — 1)  (8) 

Pl2(k/k  - 1)  = P,  2 (k  - 1/k  - 1)  + fjPl3(k  - 1/k  - 1)  + AtP22(k  - 1/k  - 1) 

+ (a  + 0At)P23(k-  1/k  - I)  + or/?P33(k  - 1/k  - 1)  + 0,2(k-  I)  (9) 

Pi3(k/k-  1)  = 7Pn(k-  1/k  - I)  + 7AtP23(k  - 1/k  - 1)  + «7P33(k  - 1/k  - 1) 

+ 0)3{k  - 1)  (10) 

P22(k/k-  1)  » P22(k-  1/k  - 1)  + 2j3P23(k-  1/k  - 1)  + d2P33(k  - 1/k  - 1) 

+ Q22(k-1)  (ID 
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Table  6.3.  Filter  Equations  for  Correlated  Measurement  Noise:  A Summary-(Continued) 


P23 (k/k)  = P23(k/k  - I)  - G(k)Kj(k)K3(k) 
P33(k/k)  = P33(k/k-  I)  - G(k)K}/k)Kj(k) 

INITIALIZATION 
x,  (0/0)  = z(0) 

*2  (0/0)  = x3  (0/0)  = 0 
P,,(0/0>  = a2(0) 

P,2(0/0)  = 0. 

P,3(0/0)  = o. 

P22(0/0)  = o|(0) 

P23(0/0)  = 0. 

P33(0/0)  = O5(0) 

CONSTANTS 
7 = exp  (-At/rM  ) 

0 s rM  ( 1 “ >) 
a = Ty  (7  + At frM  - I ) 

P = exp  (-At It) 

PROCESS  NOISE  (CONSTANT) 

Q33  = 2Ato£/TM 


(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 

■ 34) 

(35) 

(36) 


(37) 

(38) 

(39) 

(40) 


(41) 


t 
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Table  6.3.  Filter  Equations  for  Correlated  Measurement  Noise:  A Summary  (Continued) 


Q:i 

= AtQ33/2 

(42) 

Q22 

= 2*tQn/3 

(43) 

0,3 

= Q22/2 

(44) 

0,2 

= 3AtO,?/4 

(45) 

On 

= 2AtQ12/5 

(46) 

PREDICTION 

X|(t  + tp/t)=  x,(t/t)+  x2(t/t)tp  + x3(t/t)T‘(exp(~tp/rp) 

+ tp/rp-l|  (47) 

C.  RESULTS  AND  SENSITIVITY 

We  will  now  apply  this  algorithm  to  assess  the  effect  of  measurement  correlation  on 
filter  performance.  The  latest  version  of  the  filter  and  predictor,  as  developed  in  previous 
sections,  remains  the  same  ir.  form.  The  correlated  noise  filter  simply  adds  some  additional 
terms.  It  should  be  noticed  that  the  ease  of  r = 0 (white  noise)  recovers  the  identical  results 
as  before.  Figure  6. 1 presents  the  results  for  this  case  and  also  those  with  several  other 
values  of  up  to  r = 1.0  second.  As  expected,  the  presence  of  autocorrelation  in  the  measure- 
ment errors  tend  to  degrade  filter  performance  both  from  the  standpoint  of  settling  lime 
and  steady  state  values.  This  is  explained  by  the  fact  that  the  data  being  processed  is  not  as 
informative  (due  to  the  presence  of  the  random  bias)  as  independent  data.  In  fact,  the 
performance  degradation  is  rather  severe  when  the  correlation  time  is  greater  than  0.1  second 
Of  course,  the  amount  of  serial  correlation  is  a function  of  the  particular  sensor  system  and. 
therefore,  not  a parameter  under  our  direct  control.  From  a software  standpoint,  we  simply 
have  to  live  with  whatever  sensor  autocorrelation  and  subsequent  performance  degradation 
with  which  we  are  presented.  It  is  hoped,  however,  that  this  consideration  will  be  properly 
weighed,  undoubtedly  along  with  many  others,  when  decisions  as  to  a choice  of  sensor  suits 
are  effected. 

For  the  one-dimensional  work,  values  of  the  measurement  error  standard  deviation  of 
o = 5 and  20  yards  were  assumed  and  used  throughout  the  report  in  all  sections  that  deal 
with  only  one-dimensional  filtering.  In  fact,  we  have  always  normalized  our  criteria 
o,  (t  + IQ/t)  by  o as  if  o,  (t  + 10/t)  wetv  then  independent  of  a.  While  this  is  true  for 
certain  cases  (see  Section  III.  for  example),  it  is  not  exactly  true  in  general.  Therefore,  in 
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Figure  6. 1 . Correlated  Measurement  Noise  Filter  with  Various  Noise  Correlation  Times 
r Normalized  Standard  Deviation  of  Predicted  ( 10  second)  Position  Error 
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this  section  where  we  are  considering  the  effect  of  one  uncontrollable  sensor  parameter  t.  it 
seems  appropriate  to  consider  another,  i.e..  a.  in  Figure  6.2.  we  have  therefore  plotted  a few 
cases  with  different  parameters  in  order  to  sec  the  effect.  Examples  with  white  noise  (r  = 0) 
and  one  with  colored  noise  (r  = 0.30  seconds)  were  chosen.  As  can  be  seen,  the  effect  is  not 
large  but  is  not  negligible.  We  find  that  the  larger  the  value  of  o.  the  smaller  the  steady-state 
value  of  0|  (t  + 10/t).  This  indicates  that  the  actual  performance  criteria.  <J\  (t  + 10/t).  has 
less  than  a one-to-one  sensitivity  to  o. 


k - 16  Hertz 

CTjftO)2  300  yirds/seccnd 
CTj(0)=  JO  yerds/secend* 


- 0.1  y»rd$/second* 
T||  = 20  seconds  J 
Tp  - 20  seconds  I 


(a,  r) 

( (0,0.30) 
(20,0.30) 
(10.0) 
(40,0.30) 
(20.0) 

(40,0) 


TIME  f (SECONDS) 


Figure  f>.2.  Correlated  Measurement  Noise  Filter  .'.nd  Various 
Combinations  of  o and  r Normalized  Standard 
Deviation  of  Predicted  ( 10  second)  Position  Error 
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While  previously  discussing  the  question  of  a reasonable  ^ange  of  values  for  r.  we 
alluded  to  the  fact  that  we  cannot  always  expect  r to  be  the  same  constant.  In  fact,  we 
should  expect  our  estimates  of  the  sensor  statistics  to  merely  approximate  the  average 
behavior  of  the  sensor  system.  When  the  sensor  system  is  operated  under  various  conditions 
and  tuned  by  various  personnel,  we  should  expect  some  variability  from  our  estimates  to 
be  realized.  We  then  want  to  consider  the  degradation  in  filter  performance-i.e.,  the 
sensitivity  when  the  actual  sensor  parameters  are  other  than  those  assumed  by  the  filter. 

In  order  to  do  this,  we  must  again  calculate  the  value  of  the  actual  covariance  when  the  gains 
are  computed  suboptimaliy  (with  the  wrong  values  of  r and  o denoted  by  the  subscript  0 
and  compare  this  covariance  with  the  optimal  covariance.  For  the  serially  correlated  measure- 
ment error  filter,  the  actual  covariance  is  propagated  by  Equation  (2.16)  in  the  extrapola- 
tion stage  (the  same  as  the  white  noise  filtcr>  and  by  the  following  equation  for  the  update. 


PArr<k/k)  = PAtT(k/k-  1)  - |PAcr<k/k-  !)H£cT(k)-  AACT(k)|  KT(k) 
- K(k) |HAcT<k)PA{-T(k/k-  I)  * AA<-j(k)l 


(6.34) 


+ K(k)GACT(k)KT(kl 

where  AArr<k).  BA<-T(k)  and  GA<  T(k)  are  defined  in  Table  6.2  with  the  actual  values.  We 
again  use  0.  as  defined  in  Equation  (4. 1 7).  as  our  figure  of  merit  to  measure  optimality. 

In  Figure  6.3. 0 is  plotted  for  our  standard  problem  tp  = 10  seconds  integrated  for 
T = 50  seconds.  In  (a),  the  ratio  of  the  actual  to  assumed  value  of  r is  allowed  to  vary  from 
0 to  2 for  several  assumed  values  of  77 . Similarly  in  (b).  the  ratio  of  the  actual  assumed 
value  of  a is  allowed  to  vary  from  0.25  to  2.00  for  both  white  noise  (r  = 0)  and  colored 
noise  (r  = 0.30).  In  both  cases,  we  find  greater  sensitivity  when  we  overestimate  the  param- 
eters (17  > r and  0(  > o)  and  only  very  moderate  performance  degradation  when  we  under- 
estimate. In  either  case,  it  takes  a severe  error  in  the  sensor  parameter  estimates  to  seriously 
degrade  the  Kalman  filter.  This  is  somewhat  surprising  and  comforting  since  such  filters 
tend  to  rely  rather  heavily  on  a priori  estimates.  W'e  conclude  that  it  is  desirable  not  to 
overestimate  r or  a 
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Before  leaving  the  subject  of  autocorrelation  noise,  let  us  consider  one  more  question. 

It  is  interesting  to  determine  the  cost  involved  if  autocorrelated  noise  were  processed  by  a 
white  noise  filter  (rf  = 0).  With  this  information  in  hand,  we  also  want  to  consider  the 

relative  performance  of  the  effective  white  noise,  suggested  by  D’Appolito  (1971).  D'Appolito  i 

says  that  an  equivalent  white  noise  of  variance 


-'eff 


(6.35) 


“will  produce  the  same  estimation  error  as  the  original  first  order  Markov  process  of  variance 
o2 . Notice  the  factor  multiplying  o2  is  always  greater  thzn  one.”  In  Figure  6.4. 0 is  again 
plotted  for  the  case  of  rf  = 0 (white  noise  filter)  when  the  actual  noise  has  correlation  time 
r.  We  find  that,  if  the  measurement  error  standard  derivation  assumed  by  the  white  noise 
filter  o{  is  equal  to  the  nominal  value  o.  the  performance  drops  off  rather  sharply  when  * is 
greater  than  At  ( 1/16  second  here).  Unfortunately,  this  analysis  also  indicates  that  the 
effective  white  noise  (<>r  = oCfr)  performs  no  better  and.  in  fact,  considerably  worse  as  the 


Figure  6.4.  Performance  of  White  Noise  Filter  with  Various  Values  of  o 
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correlation  time  grows.  Since  increasing  the  nominal  measurement  error  does  not  help,  as 
an  experiment  the  variance  of  the  white  noise  driving  term  for  the  first-order  Markov  process 
was  tried.  We  sec  frem  Equation  (6.9).  that  om  is  less  than  the  nominal  value  o.  Thb  value 
(Of  = ow ) produces  about  the  same  results  as  for  Of  = o.  We  conclude,  therefore,  that  near- 
optimal  performance  ( 0 = 0.9).  for  the  white  noise  filter  is  obtained  only  when  r < At.  For 
r > At.  the  white  noise  filter  performance  drops  rapidly.  It  does  not  appear  from  this 
analysis  that  a white  noise  filter  with  an  equivalent  white  noise  (other  than  Of  - o)  will  help. 

A final  decision  on  the  use  of  this  correlated  measurement  noise  algorithm  must  depend 
upon  the  final  selection  of  a sensor  suite  chosen  for  the  fire  control  system. 


VII.  TECHNIQUES  FOR  THE  REDUCTION  OF  COMPUTATIONAL  BURDEN 

Several  aspects  of  the  filter  developed  to  this  point,  such  as  a desire  for  a high  data 
rate  and  rcai-timc  covariance  propagation,  begin  to  impose  a significant  computational  burden 
on  a small  digital  computer  that  might  be  utilized  in  a fire  control  system.  Further  develop- 
ment of  the  three-dimensional  filter  (in  Section  IX)  and  the  addition  of  the  dual-bandwidth 
adaptive  features  (in  Section  V')  serve  to  multiply  these  computational  requirements.  Ob- 
viously. anything  that  can  be  done  to  reduce  the  calculations  required  by  the  basic  one- 
dimensional  filter  could  l*e  significant.  In  thb  section,  two  possibilities  ate  explored. 

The  high  cycling  rate  of  the  filter  might  possibly  be  offset  by  data  compression  or 
■*prefiltering.**  U\.  processing  data  at  a higher  rate  than  the  filter  cycles  if  thb  can  be  done 
without  significantly  degrading  performance.  The  real-time  propagation  of  error  covariance 
on  a fixed-point  computer  of  limited  word  length  (say  16  bits)  poses  other  problems.  Con- 
sidering the  potential  range  of  the  elements  of  the  error  covariance  matrix,  which  involves 
the  squares  of  both  rather  large  and  rather  small  numbers  during  one  run.  it  was  found  that 
such  a computer  would  have  to  perform  the  covariance  calculations  in  double  prevision  ir, 
order  to  allow  sufficient  scaling.  There  is  also  the  possibility  of  the  loss  of  the  positive 
definite  property,  required  for  the  error  covariance,  due  ;o  numerical  error  that  results  from 
finite  word  length  calculations.  Such  a condition  b disastrous  as  it  usually  leads  to  in- 
stability and  ultimately  total  failure  of  the  filter.  The  possibility  of  this  occurring  for  our 
case  is  rather  remote,  however,  since  process  noise  is  always  added  which  tends  to  place  a 
well-defined  lower  hound  on  the  stead> -state  error  covariance.  These  problems  can  be 
eliminated  through  the  introduction  of  an  error  covariance  square  root  which  can  be  propa- 
gated in  pi.  v of  the  covariance.  There  arc  also  other  possibilities  for  determining  P.  such  as 
the  infomiaiH  •»  matrix  (the  inverse  of  error  covariance)  or  its  square  root,  but  these  will 
not  be  explored  n this  report. 

Another  ver,  important  possibility,  currently  under  investigation  and  not  included  in 
this  report,  is  the  tv*.  -Hlily  of  using  functional  approximations  for  the  error  covariance 


wn'7.  h S“Ch  Unct,0f  would  be  computed  in  the  laboratory  prior  to  implementation  and 
uid  be  functions  of  time  as  well  as  all  the  parameters  (i.e..  At,  o,  oM  . and  r„  ) The  ap- 
proach presently  oemg  studied  would  work  as  follows.  The  error  covariance  matrix ’s 
ca  culated  (via  the  Kalman  filter  equations)  as  a function  of  time  and  over  a range  of  all 
he  parameter  space  and  stored  on  the  computer.  The  steady-state  solution  is  first  fit  using 

ta  e ZT  °Ver  T™1"  range-  (W°rk  iS  a,S°  Under  to  -olv=  exactly  the  stead* 
te  nonlinear  matrix  Ricotti  equation  for  P but  a solution  has  not  yet  been  obtained.) 

m a a fZ^of  ttead>Te  S°!Uti0n  and  thC  initSal  C°nditi0nS> the  transient  Phase  * 

in  the  uZ  fit  T v reSUkmf  fUnCti°nS  are  then  USed  to  ca,cula*  ^e  gain  matrix 
m the  usual  fashion.  A sensitivity  analysis  ,s  performed  in  order  to  assure  an  accurate  ap- 

P oxlinat,on  for  the  optimal  gain.  If  this  work  is  a success,  most  of  the  computational 

ta  Z d TmlKdmm  n"Cr^hich  “ dUC  * *•  covariance  equation"  - can 

jZoZrlr*  n “ ‘*“*f  0plimisti,;  as  ‘°  the  futurc  °r  “*  method  its  appli- 
«.ation  to  the  GFCS  filtering  problem. 


A.  PRHFILTERS 


Preliltcring.  as  the  term  will  be  defined  i.i  this  report,  means  the  processing  of  data 
whirl,  is  avatlaolc  at  a rate  higher  than  that  at  which  we  wish  to  cycle  the  Kalman  filter  It 
“ . ° rr1*  r'a''TcJ  10  » *«»  compression  Suppose,  as  shown  in  Figure  7 I that  the 

“Th, ™ d t cwy  t * «« « -* ..  ^ zJt  z* 

rate  p tunes  the  fitter  cycltng  rate.  We  will  therefore  have  „ measurements,  equally  spaced 

wan,  , l’“"  mad'  S““  ‘he  laS‘  niKr  »•  Mme  ttk  - I J and  which  we 

nt  to  process  at  time  t(k).  There  are  undoubtedly  several  possible  methods  of  aggregating 

(or  lump,  ag)  or  otherwise  smoothing  these  additional  measurements.  However  w will  con 
stder  (for  now  a,  leas,;  the  simples,  effective  method  of  doing  this,  nameiy.  aver,  ginT 


es,cmfial7Z  P ° ' ^ to“"*  «*«•  "»  measurement  noise 

essentially  masks  any  time  variation  in  the  signal,  data  averaging  is  an  effective  means  of  data 

compression  w„ . small  loss  of  information.  Actually,  we  are  not  quite  (but  almoin ’the 

from"  r,  ,oP4  Hmr  1 ?' ,hiS  IT1'  T”  if  wc  wam  '»  «*«**«  measurements 

. rt /.  there  is  actually  some  floaty  information  that  could  be  extracted  from 

the  measurements.  The  variance  of  this  velocity  estimate,  however,  is  so  larges VeMvc  to 

tim,  T £.*h  ?!  ,lK' Ka' ra"  """• ,hal  "S  makes  essentially  no  improve- 

~ \ ' howevcr'  a significant  increase  in  the  computation  required  to  --rocess  such 

mcm°!7dmeaSaTT-  Wc  Wi"  tlwrefore  compuIe  “ mtmwhm,  prefiltereu  measure- 
mem  h.sed  upon  a teehmquc  similar  to  data  averaging -residual  averaging.  By  averaging  the 

Z ' : * <TSCd  t0  ,hC  "Kasura™'ts  « uccount  for  ihe  tesfima ted)  iarg* 

motion  over  the  pref liter  interval.  8 
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Measurements  available  every  At//4 


p MEASUREMENTS 


Kalman  filter  cycles  every  At 

Figure  7.  i . Prefilters  Relation  of  Measurement  and  Filtering  Timing 


The  averaged  (or  prefiltered;  residual  is  simply 


Vpftk/k-  I)  Kk  - 1 + i/Vik-  I)  (7.1) 

All  notation  in  this  section  on  prefilters  will  be  referenced  to  the  filter  cycling  rate  so  that 
times  will  be  in  fractions  (i/p)  of  the  time  between  filter  cycles.  At.  Substituting  the  defini- 
tion for  the  a priori  residual,  we  find 


np|(k,k  - I ) = — ^ (/.(k  - 1 + i/p)  - X|  (k  - 1 + i//tlk  - I )] 

^ i=i 

(7.2) 

1 i=*' 

= 7lk) ^ X|  (k  - 1 + \lp\k  ■ l) 


where 


I 1 u 

7(k)  = — V z(k  - 1 + i/p)  (7.3) 


74 


I 


is  the  average  measurement.  It  is  convenient  to  express  the  estimates  xj  at  times  (i/p)  At 
relative  to  the  extrapolated  value  at  time  t^.  This  is  done  using  the  transition  matrix  and 
can  easily  be  found  to  be 


xj  (k  - 1 + i/pik  - 1)  = Xj  (k/k  - 1)  + I j At  x2  (k/k  - 1) 
+ arj  x3  (k/k  - 1) 

a* = a[(i"  !)At] = TM  +(i_ ')  At/TM-  ] 

Substituting  Equation  (7.4)  into  (7.2)  and  taking  the  indicated  sums,  we  find 


where 


vpr(k/k-  1)  = z(k)  - Xj  (k/k  - I)  + 


At  x2 (k/k  - 1 ) 


t-n 


~hL  <*•  *3 o^/k - 1 ) 


We  therefore  define  the  effective  prefiltered  measurement  as 

zp,(k)  = z(k)  + p2  x2(k/k  - I)  - m xj(k/k-  I) 

where 


= At  (|  - _L  £ j) 
\ t*’  £1  / 


and 


**3  = r « 


i=n 


— y r /m  - i 

py 


TM  P2 


re  constants  computed  once  before  implementation.  Therefore 

t'pffk/k  - 1)  = zpr(k)  - X|  (k/k  - I) 


(7.4) 


(7.5) 


(7.6) 


(7.7) 


(7  8) 


(7.9) 


<7.i0) 
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We  note  that  if  m - I . then  ps  and  pj  vanish  and  zpf(k)  becomes  z(k)  as  before. 

Let  us  now  consider  the  effect  of  prefiltering  on  the  measurement  error  statistics.  From 
Fquation  (7.7).  it  is  obvious  that,  strictly  speaking,  an  error  in  zpf  is  a function  of  the  errors 
ei ( k/k  - l)andej(k/k-  1 ) as  well  as  the  measurement  errors.  However,  consideration  of 
the  small  values  of  pi  and  /13  reveal  these  effects  to  usually  be  negligible.  Under  this  assump- 
tion. we  find  that  the  prefiltered  measurement  error  is  simply  the  average  of  the  individual 
measurement  errors.  That  is 


ir*» 


Vpf(k)  = 


= - y v<k  - 1 + i 


i/p) 


(7.11) 


The  standard  deviation  of  the  prefiltered  measurement  error  is  (assuming  the  measurement 
error  standard  deviation  constant) 

«P,( k ) = i:  |v',(kM 


= — i: 

U' 


l-n  l = M 


v y~*  v(k  - 1 + i’p;  v<k  I + j /p) 
i=i  -1 


(7.12) 


The  term  in  the  brackets  is  the  average  autocovariance  of  all  the  measurement  error  pairs. 
I’ach  of  these  can  be  expressed  in  terms  of  the  nominal  autocorrelation  coefficient  for  lag 
At  as 


K,j  = L |V(k-  I + i/p)  V(k  - I +j  p>| 


A 


4 


= a7(k)e.\p|-  it(i>-  t(j)t/r | 


= o*(k)  exp 


= o2(k) p(k)1' 


Sub-  tituting  into  liquation  (7.12).  we  find 


(k)/o2(k) 


— L E p{k)h 

/*-  i=.  71 


(7.13) 


(7.14) 
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Notice  that,  if  p(k)  equals  zero  (white  noise),  and  noting  that  0°  is  unity.  o2r(k)/o2(k)  is 
the  familiar  1/p.  A similar  expression  can  be  written  for  the  autocorrelation  of  the  Vpf(k)’s 
which  is 


RPf(k.  k - 1)  = E lVpf(k)  Vpf(k  - 1)|  (7.15) 

Using  Equation  (7.13).  we  find 

, i=a  \-n 

RPf(k.  k - l)/o2(k)  = — £ £ p{k)  i i + "i/M  <7.16; 

Tm  hi 

In  this  case,  ifp(k)  is  zero.  RPf(k)  is  zero.  Therefore,  if  the  original  measurement  errors  are 
uncorrelated,  the  prefiltered  measurement  errors  will  also  be  uncorrelated.  Also,  if  p = 1. 
we  recover  o2,(k)  = o2(k)  = o2(k)  and  Rpr(k.  k - 1)  = o2(k) p(k). 


The  correlation  coefficient  for  the  prefiltered  measurement  errors  is 


Ppf(k)  = Rpr«k.k  l)/o2,(k) 


(7.17) 


and  the  effective  prefiltered  measurement  correlation  time  constant  can  be  calculated,  if 
desired,  as 


rpf(k)  = At/Vn  | t/ppr(k)i 


(7.18) 


In  Figure  7.2.  (he  ratios  of  oP<7o  and  rPf/r  are  plotted  as  a function  of  p (as  though  p 
were  a continuous  variable).  Since  opf/o  and  rp,/r  are  functions  only  of  p and  the  nominal 
correlation  coefficient  p.  several  values  of  t were  chosen  to  specify  p.  The  nominal  data 
rate  is  k = 16  Hertz  so  that  the  filter  cycling  time  for  each  value  of  p or,  the  graphs  is  At  = 
p/k.  We  find  that  the  ratio  opflo  is  always  < 1 for  p > 1 . ar.d  the  measurement  etror  reduc- 
tion is  greatest  for  the  white  noise  case.  The  ratio  of  rpf/r  is  > l and  increases,  on  the  other 
hand,  as  p increases.  The  ratio  decreases  as  r increases,  however. 

Let  us  now  consider  the  effect  of  prefiltering  on  filter  performance.  We  find  there  are 
three  factors  whose  interaction  we  must  consider.  The  increase  in  the  time  increment  be- 
tween filter  cycles  (a  reduction  of  the  filter  cycle  rate)  and  the  increase  in  the  autocorrela- 
tion of  the  prefiltered  error  both  tend  to  increase  covariance  and  thus  degrade  performance. 
Fortunately,  however,  the  decrease  in  the  prefiltesed  measurement  error  variance  tends  to 
decrease  the  error  covariance  and  improve  filter  performance.  The  net  effect  on  filter  per- 
formance was  found  to  be  essentially  negligible  with  less  than  2 percent  variation  in  op(t  + 
10/t)  when  r ranges  from  0 to  0.5  seconds  and  p varies  from  1 to  10.  It  has  thcrcfofe  been 
concluded  that  this  method  of  prefiltering,  based  upon  residual  averaging,  offers  us  a 
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Figure  7.2 


computationally  inexpensive  way  to  achieve  high-data-rate  filtering  with  negligible  sacrifice 
in  performance  relative  to  high-data-rate  full-Kalman  filtering. 

A small  problem  was  encountered  from  the  use  of  this  method  of  prefiltcring  which 
will  now  be  discussed  along  with  the  solution.  As  mentioned  previously,  the  residual  aver- 
aging technique  of  prefiltcring  causes  the  actual  effective  measurement  error  statistics  to  be 
a function  of  the  estimation  error,  as  can  be  observed  in  Equation  (7.7).  It  was  previously 
assumed  that  this  portion  of  the  prefiltered  measurement  error  was  negligible  as  far  as  the 
calculation  of  the  effective  prefilter  statistics  is  concerned.  This  assumption  is  valid  except 
during  the  initial  covergence  period  when  the  estimation  error  covariance  is  quite  large. 

Upon  comparing  the  calculated  covariance  with  the  actual  (simulated)  Monte  Carlo  error,  it 
was  discovered  that,  during  the  initial  covergence  interval,  the  actual  errors  consistently 
exceeded  the  calculated  error  standard  deviation  but  that  the  effect  disappeared  after  a 
couple  of  seconds.  The  author  refers  to  the  phenomenon  as  “prefilter  lag."  which  can  be  ex- 
plained by  considering  Figure  7.3.  During  the  initial  convergence  period,  the  expected 
value  of  the  residuals  evaluated  at  the  times  of  actual  measurement,  as  shown  by  the  circles, 
can  significantly  increase  over  the  prefiltering  interval  At.  The  lag  problem  occurs  because 
an  average  of  these  residuals  does  not  have  the  expected  error  level  (variance)  that  the  filter 
calculates.  The  standard  deviation  of  the  actual  residual  average  is  always  less  than  the  value 
which  the  filter  expects.  o,.(kik  - 1 ). 

After  trying  several  methods  to  rectify  this  situation,  the  author  finally  found  a tech- 
nique whereby  the  prefiltered  residual  would  be  "scaled  up"  so  that  its  variance  would 
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o - Value  Calculated  by  Error 
Covariance  Extrapolation 


match  that  calculated  by  the  filter.  This  was  done  as  follows.  First  it  was  noticed  that  the 
variation  of  a,,  over  the  interval  At  is  approximately  linear  (since  velocity  estimation  error 
dominates)  so  that  we  can  write  a:,  evaluated  at  the  measurement  points  i as 

o,.(i)  = (1  - i/p)o„(k-  l/k-  1)  (i/p)  ajk/k  - I)  (7.19) 

The  average  v;'»e  ofo,,(i)  occurs  at  the  point  i = (p  - 1 >/2  or 

o,AVt.<k)  =(~p)o,.(k-  l/k  - I)  + ^~)o,(k/k-  I)  (7.20) 

'finally,  the  factor  by  which  we  want  to  scale  the  prefiltered  residual  is  o,.(k/k  - I )/o,,AV(;(k) 

or 


PI  + r,.(k) 

_x_ 

r,.(k)  - 11 

L : 

V 

& 

where 

r,  tk)  = o,.(k  - l/k  I )'o.(k/k  - It  (7.22) 

is  the  ratio  of  residual  levels  across  the  interval.  Clearly,  if  there  is  no  prefiltering  <p  = I ) or 
if  the  filter  is  in  steady-state  operation  (r  = I >.  this  factor  is  one  and  does  not  influence 
performance.  It  was  observed,  however,  that  this  prefiltered  residua!  scaling  eliminated  the 
prefilter  !•••:  entirely  and  was  an  economical  "fix”  to  implement. 

Concerning  future  work  in  the  area  of  prefiltering,  the  author  believes  that  an  even 
greater  amount  of  data  compression  can  he  achieved  by  using  simple  data  processing  tech- 
niques such  as  least  squares.  The  autltor.  as  stated  previously,  feels  that  target  motion  is 
highly  band-limited  in  that  the  upper  limit  of  frequency  is  fairly  discernible  from  aero- 
dynamic limitations  of  aircraft  and  missiles.  Consideration  of  the  Shannon  sampling 
theorem,  with  the  usual  factor  of  ten  pul  in  to  account  for  noisy  data,  tells  us  that  it  is 
necessary  to  cycle  our  target  motion  estimator  at  only  2 to  1 Hertz.  Actually,  cycling  at 
1 Hertz  is  marginal,  according  to  Shannon,  but  we  are  probably  not  interested  in  lor  capable 
of)  actually  recovering  ail  the  very  highest  frequency  target  motion  hut  merely  want  to  track 
through  it  without  diverging.  Also,  verv-high-frequency  target  :u  celeration  usually  results 
in  very  small  actual  displacement  that  tends  to  be  unobservable  as  it  is  down  in  the  meas- 
urement noise.  A simple,  constant-velocity  least -squares  fit  over,  say.  a compression  interval 
of  1 second  and  a data  rate  of  perhaps  32  Hertz  appears  v«**y  attractive  at  this  time  and 
should  be  performance-tested  as  soon  as  possible. 
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B.  COVARIANCE  SQUARE  ROOT 


The  covariance  square  root  method  always  insures  a symmetric  positive  definite  error 
covariance  matrix.  More  significantly,  the  square  root  formulation  can  propagate  the  error 
covariance  in  single  precision  as  accurately  as  the  conventional  error  covariance  methods 
do  in  double  precision.  Basically,  this  technique  propagates  a matrix  S,  instead  of  P.  where 
S is  defined  by  the  relation 


P(k/j)  = S(  k/j  > ST  < k/j ) 

(7.23) 

This  definition  of  S is  not  unique,  however.  For  example,  in  our  situation.  Pisa(3  X 3) 
symmetrix  matrix  specified  by  6 variables.  Since  S is  also  a (3  X 3)  matrix,  for  which  9 
variables  are  required  for  specification,  we  find  ourselves  with  three  extra  degrees  of  freedom 
in  S which  might  be  used  advantageously.  We  have  chosen,  for  reasons  to  be  made  clear 
shortly,  to  complete  the  definition  by  requiring  S to  be  upper-triangular.  As  we  shall  find, 
the  definition  of  S as  upper-triangular  does  not  insure  that  S will  remain  in  this  form  when 
propagated  through  the  filter  equations.  The  upper-triangular  definition  of  S requires 

P;i  (k/j)  = Sy,(k/j)  + Sf2(k/j>  + Sf3<k/j> 

(7.24a) 

Pi2<k/jl  = Siitk/j)  Sijtk/j)  + S|j(k/j)S23(k/j) 

(7,24b) 

Pn<k/j>  = S 1 3 ( k/j  > S3  3 ( k/j  > 

(7.24c) 

p22*k/j>  = S?2«k/j)  S§3«k/j) 

»7.24d) 

P^jtk/j)  = S2,(k/j)S33(k/j) 

(7.24e) 

Pjj(k/j)  = S^tk/j) 

(7.240 

The  inverse  relations  are 

SyyCk/ji  = y/Pniklj) 

(7.25a) 

S23(k/j)  - P2 3 < k/j )/S3 3 < k/j > 

(7.25b) 

S:2fk/j)  = x/p22<k/j)  - S$3(k/j) 

(7.25c) 

S|j(k/j)  = P,3(k/j)/S33(k/j) 

(7.25d) 
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S,2(k/j>  = (P,2(k/j)  - S13(k/j)S23(k/j)l/S22(k/j) 


(7.25e) 


Sn(k/j)  = ^PnOc/j)  - Sf2(k/j)  - Sf3(k/j)  (7.250 

Now  ict  us  consider  the  filter  covariance  equations  as  given  in  Table  6.2.  The  error  co- 
variance  extrapolation  equation  can  be  written  as 

Ptk/k. - 1)  = P'(k/k  - I)  + Q(k  - !)  (7.26) 

where 

P'(k/k - I)  = <Hk.k-  I) Pfk  - 1/k  - l)*T<k,k-  I)  <7.27) 

it  is  easily  seen  that  the  extrapolation  to  P'fk/k  - I ) can  be  equivalently  accomplished  by 
the  square  root  as 


S'tk/k  - I)  = # k.k-  l)S(k-  1/k  - 1) 


(7.28) 


Applying  the  transition  matrix  lor  our  case,  we  find  that,  if  S(k  - 1/k  - I ) is  upper- 
triangular.  then  Sf(k/k  - I ) is  also  upper-triangular  with  the  following  values. 

S',t Ik/k  - 1)  = S,,(k-  1/k  - I)  (7.29a) 

s;:(k/k-  I)  = Sl2(k-  1/k-  I)  + AtS22(k-  1/k  - I)  (7.29b) 


S’,3<k/k-  I)  = S,,(k-  1/k  - I)  + At  S23(k  1/k  - I) 


+ aSjjdc  1/k  - I) 


(7.29c) 


S22(k/k  - I)  = S22(k-  1/k  - 1)  (7.29d) 

S2;,(k/k  - I)  = S:3(k-  1/k  - 1)  + 0S33(k-  1/k  - I)  (7.29e) 


S33 « k/k  - l > = y Sjjtk-  1/k  - 1) 


(7.290 


If  process  noise.  Q.  is  present,  the  addition  of  this  term  poses  a more  difficult  problem 
in  the  square  root  covariance  formulation.  There  are  several  possible  methods,  both  exact 
and  approximate,  discussed  in  the  literature.  The  simplest  and  usually  the  fastest  of  the 
exact  methods  is  the  so-called  root-sum-square  (RSS)  operation  discussed  by  Carlson  ( I973i 
an  ;>ng  others.  This  technique  essentially  recalculates  the  error  covanance.  using  Equation 
(7.24).  adds  the  process  noise  and  then  takes  the  matrix  square  root  to  determine 
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S(k/k - i)  = (S'(k/k  - DS'TOc/k-  l)  + Q(k-  1)11/2 


(7.30) 


This  last  operation,  using  Equation  (7.25).  is  commonly  known  as  Choiesky  decomposition 
in  triangular  form.  An  alternative  exact  method  is  the  Householder  triangularization  algo- 
rithm. discussed  by  Kaminski.  Bryson  and  Schmidt  (1971)  which  maintains  double  preci- 
sion accuracy  but  with  higher  computational  cost.  Schmidt  (1970)  also  discusses- but  does 
not  recommend-  a technique  which  requires  two  matrix  inversions  and  results  in  a loss  of 
the  upper  triangular  form.  In  addition  to  considering  these  exact  algorithms,  the  author 
also  found  the  technique  of  Wu  (1973)  very  interesting.  Wu  assumed  a square  root  addi- 
tive process  noise  of  a form 


S(k/k  - 1)  = S'(k/k  - l)  + q'(k-  1) 


(7.31) 


where  we  might  think  of  q'(k  - I ) as  equivalent  in  some  sense  to  Q,  /2(k  - I).  Equation 
(7.31)  is  admittedly  not  equivalent  to  the  original  extrapolation.  Equation  (7.26).  since  two 
additional  terms  appear  in  the  error  covariance.  That  is 


P(k/k  - 1)  = S'(k/k  - 1)  S^fk/k  - 1)  + q'(k  - Dq^tk-  I) 
+ S'fk/k  - I ) q'T(k  - 1)  + q'lk-  DSTfk*-  I) 


(7.32) 


Wu's  principal  argument  for  this  assumption  is  that  the  process  noise  itself  (Q>  is 
rarely  known  exactly  in  the  first  place,  being  “basically  empirically  determined  data.”  The 
method  is  also  no  different,  in  principle,  from  the  epsilon  technique  in  which  a somewhat 
arbitrary  process  noise  is  added  to  prevent  divergence.  Since  divergence  prevention  and 
control  of  the  filter  bandwidth  are  the  primary  uses  of  process  noise  in  our  application, 
this  argument  does  not  seem  unreasonable.  The  main  difficulty  in  Wu's  method  is  in  appli- 
cation. It  is  very  difficult  to  control  the  amount  of  process  noise  which  we  would  like  to 
add.  Notice  that  the  last  two  terms  involve  S'  so  that  the  amount  of  process  noise  actually 
added  is  a function  of  the  current  error  covariance.  These  two  terms  are  definitely  not 
negligible  since  each  is  usually  larger  than  the  q'  q'T  term.  The  reason  for  this  is  that  S'  is 
usually  much  larger  than  q\  The  author  found  that,  by  trial  and  error,  it  is  possible  to  find 
a q'  (much  less  than  Q,/2 ) that  produces  essentially  the  same  steady-state  error  covariance 
;»s  the  covariance  filter  for  a given  situation.  Unfortunately,  no  relation  can  be  found  be- 
tween q'  and  Q w.iich  we  might  use  to  maintain  adequate  control  over  the  filter  bandwidth. 
Wu’s  method  was  therefore  abandoned  in  favor  of  an  exact  method. 


it  was  therefore  decided  to  use  the  RSS  technique  which  is  She  most  efficient  of  the 
exact  methods.  The  equations  for  adding  process  noise  are  then  (letting  W = Sfk/k  - I ) for 
notational  convenience) 


f 
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4 


! 


W33 

= v/s&  + 

Qj3 

(7.33a) 

W23 

= (S'23  s'33 

+ Q23VW33 

(7.33b) 

W>2 

->/§ T? 

S&  022  - Wjj 

(7.33c) 

W,3 

= (S',  3 S33 

+ Ot3)/W33 

(7.33d) 

W,2 

= (S',  2 S'22 

+ S',  3 S'23  + Q,2  - W,3  W23/W22 

(7.33e) 

w,, 

= v^S',2,  + 

s;22  + sft  + On  - W|2  - w?3 

(7.330 

Actually,  it  might  appear~at  least  it  did  to  the  author- that  the  requirement  to  calculate 
the  covariance  in  the  RSS  method  would  undermine  the  expressed  desire  to  retain  double- 
precision  accuracy.  The  accuracy  is  not  lost,  nowever.  since  in  the  type  of  fixed-point 
computer  in  which  this  is  implemented,  a double-length  word  results  automatically  when- 
ever two  single  length  words  arc  multiplied  and  a double  length  word  must  be  used  in  order 
to  obtain  a single-length  word  from  the  square  root.  The  net  result  is  that,  while  double 
length  words  appear  in  the  intermediate  calculations,  all  the  calculations  are  really  in  single 
precision  and  very  efficient. 

Let  us  now  consider  the  update  equation.  The  original  technique  for  measurement 
update  of  the  white  noise  filter  was  developed  by  Potter  (1963)  for  the  case  of  scalar  meas- 
urements. Unfortunately,  the  Potter  technique  is  not  applicable  to  the  form  of  the  co- 
variance  update  for  the  correlated  measurement  noise  filter.  The  Potter  form  also  results  in 
a loss  of  triangularity  even  for  the  white  noise  case.  We  will  therefore  again  utilize  the  RSS 
method  suggested  by  Carlson  in  order  to  maintain  double-precision  accuracy.  The  co- 
variance  update  equation,  repeated  from  Tabic  6.2,  is 

P(k/k)  = WWT  - K*(k>G(k)  K*T(k) 

Applying  the  modified  RSS  technique  to  our  equations,  we  find  (again  letting  S = S(k/k)  for 
convenience) 


S33  'Al,  - GKJ2 

(7.34a) 

S23  = (W2J  W,:  - G KJ  K* )/S33 

(7.34b) 

S22  = + W2j  - GK?2  - 

S23 

(7.34c) 
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<7.34d) 


Si 3 = (W,3  W33  - GKf  KJ)/S33 
S;:  = (W 1 2 W22  + W,3  W23  - GKfKj  - Si 3 S23)/S22  (7.34c) 

Sm  = >/Wll  + W12  + W!3  - CK*2  - s?2  - S?3  (7.340 

Hquations  (7.33)  and  (7.34)  will  now  replace  the  covariance  equations  (8)  through  ( 13)  and 
(22)  through  (27)  of  Table  6.3.  The  gains  are  calculated  in  the  same  manner  except  that 
Liquation  (7.24)  must  be  used  to  calculate  the  required  covariances.  The  expressions  for 
the  A's  are  somewhat  simpler  in  square  root  form,  i.e., 


y\  1 ( k ) - 
+ 

A2(k)  = 
A 3 ( k ) = 


p(k)  [Sf,  (k/k  - 1)  + S j 2 (k/k  - l)Si2(k-  1/k-  1) 


(7.35a) 


S',  3 (k/k  - l)S13(k-  1/k  - 1)1 

p(k)  [S|2(k  - 1/k-  t ) S22(k/k  - 1)  + S,3(k-  1/k-  l)S23(k/k-  1)1 

(7.35b) 

p(k)  Sj  3 < k - 1/k  - 1 ) S33 (k/k  - 1)  (7.35c) 


Initialization  of  the  square  root  covariance  is  accomplished  by  anothei  USS,  similar  to 
Equation  (7.25). 

In  summary,  we  have  found  it  is  possible  to  reduce  the  computational  burden  by  a 
factor  of  p by  prefiltering  without  loss  of  performance.  Typically,  we  might  choose  ju  = 4 
and  thereby  process  data  at  a rate  of  16  Hertz  but  only  cycle  the  filter  at  4 Hertz.  By  using 
the  square-root  covariance  filter,  it  was  possible  to  eliminate  double-precision  covariance 
calculation,  thereby  significantly  reducing  the  time  required  for  covariance  computation. 

In  fact,  a comparison  was  performed  for  the  Navy  standard  mini-computer,  the  16-bit 
A/N  UYK  20,  to  estimate  the  difference  between  the  error  covariance  equations  imple- 
mented in  double  precision  and  the  square-root  covariance  equations  in  single  precision. 

The  difference  was  more  dramatic  than  anticipated  as  the  square-root  covariance  resulted 
in  a reduction  of  computer  time  by  a factor  of  4.56!  This  will  translate  into  an  overall  filter 
computation  reduction  of  a factor  slightly  less  than  four-since  the  filter  obviously  does 
other  things  besides  covariance  calculations.  The  combination  of  both  these  techniques 
reduced  the  required  filter  computation  time  by  a factor  of  12  to  16  over  the  high-data-rate 
covariance  filter.  We  find  that  these  reductions  are  well  advised  since  the  dual-bandwidth 
adaptive  filter  introduces  an  increase  factor  of  two  and  the  three-dimensional  filter  (Section 
IX)  an  increase  factor  of  three.  That  is,  we  will  be  (approximately)  able  to  run  a three- 
dimensional,  double  adaptive  filter  with  prefiltering  and  square  root  covariance  in  one  half 
the  time  of  single  one-dimensional  error-covariance  filter  without  prefiltering. 
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VIH.  '-.^DICTION 


This  section  concerns  itself  with  th;  qwt  tion  of  the  prediction  of  future  target  posi- 
tion based  or:  current  estimates  of  position  -nc  rates.  We  will,  for  the  moment,  divorce 
ourselves  from  the  target  modeb  utilized  fot  the.  filter  and  reconsider  modeb  foi  the  long- 
time prediction  problem  (i.e.,  tp  = 5 to  30  seconds).  We  will  think  of  the  filter  only  as  a 
"black-box"  algorithm  which  provides  as  outputs  estimates  of  current  position,  velocity, 
acceleration,  and  associated  error  covariance.  Thb  information  will  serve  as  input  to  the 
predictor.  The  question  b how  the  predictor  can  best  use  thb  information.  The  predictor, 
of  course,  is  the  critical  factor  in  a gunfire  control  system  in  thnt  the  calculated  future  target 
position  b the  aim  point  for  the  gun. 

The  conventional  method  of  prediction  b to  simply  extrapolate  the  assumed  dynamics 
model  used  in  the  filter  itself.  For  example,  a second-order  derivative  polynomial  filter 
(with  appropriate  divergence  prevention  modifications)  that  provides  estimates  of  the 
(possibly  varying)  acceleration  would  predict  the  position  at  time  tp  in  the  future  as 

xA(t  + tp/t)  = Xj  (t/t)  + x2(t/t)tp  + xj(t/t)  tp/2  (8.1) 

where  x (t/t)  is  the  current  state  vector  estimate.  It  lias  been  suggested  by  some  people  that 
perhaps  ihe  estimated  acceleration  should  not  be  used  in  the  prediction  even  though  such 
an  estimate  might  be  available  from  the  filter.  Thb  argument  b based  on  the  indisputable 
contention  that  ore  cannot  expect  a target  to  maintain  a constant  acceleration  over  long 
prediction  time  tntervab.  Obviously  if  the  acceleration  is  ignored  in  the  predictor,  we  simply 
predict  tangent  to  the  curve  and  have  the  constant  velocity  predictor. 

M*  + tp/t)  = x,  (t/t)  + x2  (t/t)  t„  (8.2) 

Thb  predictor  is  currently  more  prevalent  in  operational  gunfire  control  systems. 

Thb  technique,  of  course,  essentially  surrenders  any  hope  of.  at  ieast  partially,  dealing 
with  an  accelerating  target.  One  is  therefore  faced  with  a choice  of  predictors,  and  it  must 
be  conceded  at  thb  point  that,  in  the  general  case,  one  has  no  way  of  knowing  what  the 
target  will  do  in  the  future.  It  might  very  well  follow  either  Equation  (8.1 ) or  (8.2)  or  (more 
likely)  neither.  The  prediction  time  b an  important  consideration  here.  For  short  prediction 
times  (say  i to  2 seconds!,  it  docs  not  make  a great  deal  of  difference  since  very  few  targets 
can  maneuver  at  a level  to  seriously  alter  their  course  in  such  a short  time.  For  longer  pre- 
diction times  (say  2 to  20  seconds),  typifying  the  possible  air  target  projectile  flight  times 
observed  in  large  caliber  gun  systems,  it  b well  known  that  an  air  target  can  alter  its  course 
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considerably.  Indeed,  it  is  the  author's  contention  that,  in  an  engagement,  one  should  ex- 
pect the  target  not  to  maintain  a constant  velocity  or  acceleration  over  such  durations. 

The  above  reasoning  led  rather  naturally  to  a consideration  of  methods  of  somehow 
weighting  the  acceleration  in  the  prediction  equation.  The  weighting  parameters  required 
could  conceivably  be  estimated  by  considering  the  types  of  target  trajectories  in  a system 
scenario  to  determine  parametric  values  that  could  improve  prediction  accuracy.  One  ap- 
proach to  this  problem  can  be  effected  if  we  again  consider  the  target’s  acceleration  to  be 
modeled  as  an  exponentially  autocorreiated  random  variable  (as  in  Section  IV)  and  therefore 
of  a nondeterministic  nature.  Such  a model  would  conceivably  contain  any  target  maneuver, 
be  it  the  exercise  of  deterministic  strategies,  evasive  maneuvers,  acceleration  perturbations 
due  to  winds  and  turbulence,  etc.,  as  simply  a sample  trajectory  as  long  as  the  trajectory 
acceleration  profile  can  reasonably  be  described  by  the  chosen  statistics. 

Let  us  first  consider  the  prediction  transition  matrix.  If  we  again  assume  that  the 
actual  state  obeys  the  relation 

x(t  + tp ) = ‘Ktp)x(t)  + w(t  + tp.t)  (8.3) 

then,  if  we  are  given  an  estimate  x(t).  the  expected  value  of  x(t  + tp)  is  simply 

x(t  + tp/t)  = 4>(tp)x(t/t)  (8.4) 

since  w is  a zero  mean  random  variable.  Notice  that  we  have  used  additional  notation  on  w 
to  indicate  that  it  is  integrated  over  tp . It  is  quite  easy  to  confirm  that  the  transition  matrix 
also  yields  the  minimum  variance  linear  prediction.  Suppose  we  assume  some  transition 
matrix  (say  A)  which  we  want  to  determine  so  as  to  yield  the  minimum  prediction  error 
variance,  i.e..  we  assume 


x(t  + tp/t)  = A(tp)x(t/t) 

The  prediction  error  is  then 

e(t  + tpit)  = x(t  + tp|t)  - x(t  + tp) 

= A(tp ) x(t/t)  - <f>x(t)  - w(t  + tp.f.) 
and  the  prediction  error  covariance  is 


(8.5) 


(8.6) 
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P(t  + tpit)  = E |e(t  + tp It)  eT(t  + tp|t)| 

= (A(tp)  - <l>(tp)l  ft(0)  [ A(tp ) - 4>(tp)lT  (8.7 

+ A(tp)  P(tit)  AT(tp)  + Q(tp) 

The  matrix  S2(0)  is  the  covariance  of  x where 

J2(T)  = E |x(t)  xT  (t  + T)j 

is  the  autocovariance  of  x which  can  be  related  to  the  transition  and  process  noise  matrices 
by  taking  moments  of  Equation  (8.3),  i.e., 

fi(0>  - 4>(t„)n(0)4>T(tp)  = CHt„)  (8.8 


i2(tr)  = «h(tp)n(0)  (8.9] 

Using  the  techniques  described  in  Gelb  (1974).  we  construct  the  cost  function 

J = E leT (t  + tp/t)  S €(t  + tp/t)J  (8.10! 

where  S is  any  positive  semidefinite  matrix.  Choosing  S = 1,  we  find 

J = trace  | P(t  + tp/t )1  (8.11: 

It  is  a relatively  simple  exercise  to  confirm  that  the  necessary  and  sufficient  condition  to 
minimize  J is  to  choose 

A(tp)  = ♦(tp ) (8.12’ 

We  find,  therefore,  that  the  minimum  variance  predictor  is  also  given  by  the  transition 
matrix.  The  position  prediction  equation  is  then 

x,  (t  + tp/t>  = x,  (t/t)  + x2(t/t)tp 

(8.13 

+ x3(t/t)rp  lexp (-  tp/fp)  + tp/rp  - II 
which  is  the  double  integral  of  the  prediction  a<  ..deration 


xjd  + tp/t)  = x3  (t/t)  exp  (-  tp/tp) 


Notice  that  we  have  introduced  another  parameter,  rp,  which  is  the  acceleration  time  con- 
stant we  want  to  use  for  prediction.  This  value  does  net  necessarily  correspond  to  any  pre- 
vious parameters  since,  as  mentioned  in  previous  sections,  the  filter  does  not  assume  any 
particular  value  of  the  maneuver  frequency  but  instead  assumes  an  upper  and  lower  bound 
of  a range  of  values.  The  author,  at  one  time,  attcir.^l  d to  determine  a value  of  rp  on  line 
(i.e.,  during  real-time  execution)  or  adaptively,  but  to  i.-ate  no  satisfactory  technique  has 
been  found.  It  still  appears  feasible,  however,  to  estimate  an  appropriate  value  of  rp  based 
upon  the  state  of  convergence  of  the  adaptive  filter  and  additional  work  in  this  area  might 
be  productive. 

In  the  absence  of  such  a technique,  however,  it  appeared  more  promising  to  use  an 
expected  value 


rp  = E [r<*  | (8.15) 

where  the  expectation  is  taken  over  the  target  scenario.  At  this  time,  a value  of  rp  = 5 
seconds  appears  most  reasonable  to  the  author.  It  should  be  noted  that  the  actual  value  o,' 
the  predictor  is  rather  sensitive  to  the  choice  of  rp , as  can  be  seen  in  Figure  8. 1 . The  ac- 
celeration at  t = 15  seconds  is  assumed  to  be  known  as  4 yards/second2 . By  choosing  dif- 
ferent values  of  rp  (from  zero  to  infinity),  we  find  that  the  actual  predicted  position  (say 
for  t?  = 15  seconds)  varies  by  several  hundred  yards.  If  rp  is  zero,  we  essentially  are  using 
constant-velocity  prediction,  and  conversely,  if  rp  becomes  very  iarge  (or  frequency  ap- 
proaches zero),  we  predict  a constant  acceleration  path.  For  values  in  between  zero  and 
infinity,  the  predicted  acceleration  decays  exponentially,  and  th-  •’ath  eventually  approaches 
a constant  velocity - but  not  the  same  path  as  if  rp  = 0.  Notice  that  the  exponential  accelera- 
tion equation  always  predicts  a path  somewhere  between  constant  velocity  and  constant 
acceleration  and  can  never  cross  these  bounds. 

A study  conducted  by  the  author.  Clark  (1973).  compared  the  three  predictors 
(Equations  8.1, 8.2,  and  8.13)  discussed  in  this  section.  The  results  indicated  that,  while 
for  any  given  trajectory,  either  constant  velocity  or  constant  acceleration  might  be  the 
best,  the  exponential  acceleration  predictor  was  consistently  the  ove.o’i  best  (considering 
root-sum-square  prediction  error).  This  result  should  be  expected  since  the  exponential 
acceleration  matches  the  acceleration  autocovariance  of  the  scenario  much  better  than  the 
other  two  predictors. 

Before  continuing  the  filter  development,  it  is  intcresting-and  very  disheartening  to 
consider  realistically  the  type  of  prediction  error  likely  to  be  experienced  by  a gunfire 
control  system  against  a "randomly"  maneuvering  target  under  conditions  of  perfect  filter- 
ing (where  we  assume  the  current  filter  error  covariance  vanishes).  The  standard  deviation 
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Figure  8. 1 . Random  Acceleration  (Exponential)  Prediction:  Examples 


of  the  predicted  position  error  is  due  only  to  extrapolated  process  noise.  That  is 

0|(t  + tp/t  )=vT5n7  (8.16) 

where  Qi  |P  is  calculated  from  Equation  (4.16).  In  Figure  8.2.  this  quantity  is  plotted  as  a 
function  of  prediction  time  tp  for  a low.  medium,  and  high  set  of  the  parameters  om  and 
r\| . These  results  are  consistent  with  those  obtained  hy  simulation  against  targets  with  the 
same  maneuvering  characteristics.  When  we  consider  these  values  lit  li^M  of  the  ctfecrive 
radii  of  lethality  of  various  projectile/target  combinations,  we  must  conclude  that  the  ef- 
fectiveness of  such  gun  systems  (using  unguided  projectiles)  at  long  ranges  is  certainly 
questionable.  In  fact,  this  unfortunate  situation  is  the  primary  factor  limiting  the  effective- 
ness of  current  conventional  gunfire  control  systems. 

The  eventual  deployment  of  a terminally  guided  projectile  with  significantly  larger 
acquisition  "baskets”  promises  a corresponding  increase  in  the  effectiveness  of  gun  systems 
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against  maneuvering  targets.  The  random  acceleration  model,  as  developed  in  this  paper,  is 
not  intended  to  be  a “solution”  to  the  gunfire  control  prediction  problem.  It  is  really  only 
a more  realistic  approach  to  the  problem,  as  opposed  to  the  derivative  polynomial  models, 
and  it  is  felt  an  improvement.  The  author  and  his  associates  are  currently  investigating  new 
prediction  techniques  that  hold  the  promise  of  significantly  improved  prediction  accuracies 
for  particular  targets.  A report  on  these  techniques  is  expected  to  be  published  in  the  near 
future. 


IX.  THREE-DIMENSIONAL  CARTESIAN  FILTERING 

All  discussion  to  this  point  on  filtering  and  prediction  have  dealt  with  one-dimensional 
estimation  of  target  motion.  The  rationale  for  this  approach  was  to  maintain  a level  of 
simplicity  in  the  presentation  as  long  as  possible.  Obviously,  the  multidimensional  aspects 
of  the  problem  actually  had  to  be  considered  at  all  times  during  the  development.  In  this 
section,  we  will  concentrate  our  attention  on  this  subject. 

Several  factors  relating  to  the  application  must  be  considered  in  selecting  a coordinate 
system  for  (Uterine  and  prediction.  Primary  consideration  must  be  given  to  accuracy,  since 
this  is  a ilre  control  problem,  and  computational  simplicity,  since  the  algorithm  must  be 
implemented  within  the  constraints  of  a given  computer.  The  sensitivity  to  various  non- 
linear effects,  we  shall  find,  drives  the  selection  based  upon  these  factors.  Special  types  of 
measurements  (such  as  Doppler  range  rale)  or  the  existence  of  multiple  rates  between  the 
polar  measurements  would  influ  nee  this  decision  but  are  not  present  for  this  application. 

A factor  which  is  important  in  this  instance  is  the  reference  coordinate  system  required  for 
actual  filter  input  and  predictor  output.  The  requirement  for  stabilized  (pitch  and  roll 
corrected)  measurements  with  a common  reference  point  on  output  can  be  effected  with 
fewer  and  simpler  transformations  in  the  Cartesian  frame,  but  this  factor  is  not  of  primary 
interest.  The  first  two  factors,  accuracy  and  computational  burden,  were  therefore  the 
driving  forces  for  this  study  and  will  therefore  be-  considered  in  depth  in  this  section. 

The  study  presented  in  this  section  is  based  on  two  premises: 

i I > Measurements  of  target  position  are  obtained  in  spherical  polar  coordinates,  i.e.. 
slant  range  (R>,  bearing  angle  from  north  <B)  and  elevation  angle  %¥).  it  is  assumed  that,  as 
designer  of  a filier/predietor  for  a given  fire  control  system,  we  have  avaibble  “reasonably 
good”  estimates  of  the  quality  of  these  measurements,  i.e..  measurement  error  variances 
and  correlations  in  the  original  (polar)  measurement  frame.  It  is  also  assumed  that  the  de- 
signer krows  the  availability  of  this  data  such  as  the  range  of  possible  data  rates,  the  volu- 
metric coverage  of  the  sensor,  any  spatial  variation  of  the  measurement  error  statistics  (due 
to  such  effects  as  multipath),  and  estimates  of  any  target  induced  errors  such  as  glint  and 
scintillation. 


The  polar-Cartesian  transformation  equations  are: 


_x“ 

~r  sin  B cos  E” 

4c  = 

y 

= 

r cos  B cos  E 

_z  _ 

j sin  E 

~ Sc  ir ) 


(9.1) 


where  x is  east,  y is  north  and  z is  vertical.  The  corresponding  inverse  relation  is: 


X 

yjxl  + y?  + zi 

4 1 ~ 

B 

= 

arc  tan  (x/y) 

_E_ 

..arc  sin  (z/r)  _ 

= Zr(xc) 


(9.2) 


(2)  The  target  is  assumed  to  be  modeled  "best"  in  Cartesian  coordinates.  By  this,  we 
mean  that  the  target  is  more  closely  linear  and  well  behaved  in  Cartesian  <-oordinates  than 
in  polar  coordinates.  For  example,  let  us  consider  simple  linear  target  motion  which  is 
canonical  in  Cartesian  coordinates. 


dx/dt  = vx 

= Constant 

(9.3a) 

dy/dt  = vy 

- Constan* 

(9.3b) 

z = 0 

* Constant 

(9.3c) 

Using  liquation  (9.2).  we  find  that  second  (and  all  higher)  derivatives  appear  in  the 
polar  frame  as 


d2  r/'dt2 

* r (dB/dt)2 

(9.4a) 

d2  B/dt2 

= - 2 (dB/dtXdr/dt)/r 

(9.4b) 

dB/dt 

* (x  ' vx  - y • vy)/f2 

(9.4c) 

The  author  refers  to  these  second  derivatives  as  "apparent  accelerations"  and  they  ar  * 
described  in  the  literature,  such  as  Monzingo  (1972)  who  calls  them  "psuedo-maneuvers" 
and  Cantrell  (1973).  These  accelerations,  if  viewed  in  the  polar  frame.  must  either  be 
modeled  and  propagated  nonlineariy  by  the  filter/predictor  or.  even  worse,  tracked  adap- 
tively. Notice  aho,  in  Equation  (9.4).  that  any  crossing  component  can  produce  very  high 
angular  rates  if  r becomes  small.  There  arc,  of  course,  corresponding  polar  canonical  targets 
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such  as  motion  along  a ray  or  circular  motion  centered  at  the  origin.  Examination  of  the 
target  scenario  indicates,  however,  that  such  motion  would  be  encountered  much  less  fre- 
quently than  the  nominal  Cartesian  motion. 

Perhaps  the  “best”  target  coordinate  system  of  all  would  be  one  translating  with  the 
target  and  oriented  appropriately  along  the  (changing)  velocity  vector  and  two  normal  ac- 
celeration directions.  Such  a coordinate  system  would  be  very  cumbersome  to  implement, 
however,  and  would  not  change  the  fact  that  a basic  nonlinearity  exists  between  the  meas- 
urement frame  and  the  target  motion  frame.  We  will  therefore  assume  that  the  Cartesian 
random  acceleration  target  model,  where  the  maneuver  level  and  frequency  lie  within  the 
previously  described  bounds,  represents  the  true  target  and  proceed  on  that  basis.  To  form 
the  Cartesian  target  model,  an  additional  assumption  of  independence  of  target  maneuver 
between  channels  was  made.  That  is.  we  define  the  three-dimensional  Cartesian  state  vector 

* ixyz)r  <9-5) 

where  x.  v.  and  z are  each  onc-dimcnsional  three-clement  state  vectors  (position,  velocity,  and 
acceleration)  governed  by  the  radom  3cceleratK  n model  of  Section  IV.  The  transition 
matrix  is  then 


>00 

♦jik  * | 0 $ 0 

[0  0 p 

where  0 is  defined  by  Equation  (4.8)  and  the  process  noise  is 

'0  0 0 


OjDt  “ 


0 Q 0 

0 0 0. 


19,7) 


where  Q is  given  by  Equation  (4. Li).  The  assumption  of  maneuver  coordinate  independence 
Is  weak  but  depends  upv the  parri  :ufar  type  of  target  and  its  particular  angular  orientation 
in  the  coordinate  system.  The  problem  can  probably  be  resolved  by  modeling  the  target  in 
the  moving  target  freme  (mentioned  previously)  or  by  cross-correlation  studies.  Clearly, 
the  target  motion  ana'ysis  problem  is  one  that  requires  much  additknval  attention. 


Starting  with  premises  < I ) and  (2).  it  becomes  immediately  obvious  that  there  are  a 
startling  number  of  different  Hirer  configurations,  both  optimal  and  subopiimal.  that  one 
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can  construct.  Primary  classification  of  these  configurations  will  be  by  coordinate  system  of 
the  state  vector.  We  then  considered  three  classes:  (A)  Cartesian:  (B)  spherical  polar:  and 
(C)  hybrid  (a  combination  of  both).  The  details  will  be  presented  here  only  for  the  Cartesian 
portion  of  the  study  since  a suboptima!  version  of  this  class  was  ultimately  chosen.  The 
performance  criterion  was  again  chosen  to  be  the  average  (integrated)  predicted  position 
optimality  ratio 


let  »3D  Acnwl(tk  +tp/tk) 


(9.8) 


where  ojp  is  the  root-sum-squarc  of  the  three  predicted  position  error  coordinate  values. 
Based  on  the  premises  of  this  study,  we  will  now  proceed  to  form  filter  configurations  and 
test  the  performance  of  each. 


A.  OPTIMAL  NONLINEAR  FILTER 


Given  the  Cartesian  target  model  which  is  obviously  linear  in  the  Cartesian  frame,  the 
use  of  polar  measurements  implies  that  the  measurement  update  step  will  contain  the  non- 
linearity. Now  we  will  model  this  nonlinear  measurement  and  construct  an  optimal  non- 
linear filter.  The  three-dimensional  polar  measurement  is 


Zp  * |r  B E|r  + v 


K vB  vt|T 


with  measurement  error  variance 


Rr  = 


0 

0 


0 

0 


(9.9) 


(9.10) 


assuming  no  cross  correlation  of  polar  errors.  Equation  (9.9)  is  clearly  of  the  general  non- 
linear form 
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z(k)  - h(x)  + V(k) 


(9.11) 


which  can  be  treated  by  using  the  "iterated  extended  Kalman  filter.'’  The  extended  Kalman 
filter  linearizes  h(&)  about  the  extrapolated  estimate  and  uses  the  Jacobian  to  obtain  the 
filter  update.  The  iterated  filter  simply  repeats  the  process  each  time,  linearizing  about  the 
most  recent  estimate.  Details  of  this  general  technique  will  not  be  presented  in  this  report 
but  can  be  found  in  Jazwinski  (1970).  An  important  point  is  that  the  iterated  extended 
Kalman  filter,  if  nendivergent.  tends  to  remove  those  estimation  errors  which  arc  caused 
by  systematic  and/or  observation  nonlincarities.  By  iterating  until  the  changes  between 
successive  estimates  and  covariances  become  arbitrarily  small,  one  approaches,  as  closely  as 
desirable,  the  optimal  estimator  for  the  nonlinear  situation.  This  is  the  technique  utilized 
in  this  study  to  establish  the  optimal  or  standard  against  which  we  will  compare  various 
admittedly  suboptima!  (liters. 


1).  COUPLED  LINEARIZED  ITLl  HR 

Another  approach  for  three-dimensional  Cartesian  filtering  involves  utilizing  the 
polar.'CartcMan  transformation  to  nonlineariv  combine  the  polar  observations  to  form 
Carte-law  measurements.  The  three-dimensional  measurements  are  given  as 


~ r sin  B cos  E ~ 

!s  = 

y 

= 

r cos  B cos  E 

- r sin  E 

(9.12) 


with  observation  matrix 


Hr  = 


I 0000000  0 


000  I 00000 
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(9.13) 


A problem  arises  since,  although  we  input  the  exact  nonlinear  measurement,  it  is  dif- 
ficult or  impossible  to  compute  the  exact  nonlinear  measurement  error  variance.  !n  order  to 
estimate  R corresponding  to  Equation  (9.12).  we  resort  to  linearizing  the  measurement  error. 
The  sensitivity  to  the  approximation  error  of  R should  not  be  high  as  we  found  in  the  meas- 
urement error  scnsitivi'.y  ..u-^es  in  figure  6.3(b).  To  form  this  linear  approximation,  we 
first  must  linearize  th.  ■ rjtv.urcment  error.  Differentiating  Equation  (9.1 2)  and  replacing 
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the  differentials  by  "deltas”  representing  the  finite  observation  errors,  we  find  that  the 
linearized  Cartesian  errors  are  equivalent  to  a simple  rotation  of  the  polar  errors. 


where  the  Jacobian  is 


3d. 


= T (x)  Vp 


"sBcE 
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0 

0 _ 

(9.141 
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using  s and  c to  represent  sin  and  cosine. 

Again  assuming  no  cross  correlation  of  poiar  errors,  the  linearized  Cartesian  measure- 
ment error  variance  is 


Res.  = TU)  Rp  T1  (x)  (9.|f»> 

For  reference  purposes,  the  elements  of  R<-t  will  be  written  out. 

a * = sB2cE 2 o2  + r2cB2  of  + r2sB2sE2  Oj2-  (9.1?3) 

oj  = cB2cF2  of  + r2sB2  of  + r2cB2sF2  of  (9.1 7b) 

aj  = sF.2  of  +■  r2cH2  of.  (9.!7e> 


ox>  = Oy\  ~ sBtfMof  cF2  - r2  d|  + r2sF2  o,2  ) (9.1 7d) 

0\t  ~ o*\  ~ sBstvE  iof  - r2  of ) (9.1 7c) 

oy,  - o,y  = cBsF.cH  iof  - r2  op  I S9.PD 


It  should  be  noted  that  there  is  a slight  discrepancy  between  Equations  (9. 1 6)  and  (9.17). 
Every where  in  Equation  (9. 1 7)  that  of  appears,  there  originally  was  of  cos2  t as  calculated 


97 


from  Equation  (9. 16).  The  reason  the  cos2  E was  omitted  h that  tile  quantity  o\  cos2  E 
actually  is  the  value  that  remains  constant  in  the  tracking  system  over  the  range  of  elevation 
angles  from  0 to  90  degrees.  The  angular  resolution  is  really  constant  in  the  dish  frame, 
i t.,  in  the  polar  coordinate  system  that  is  aligned  with  the  sensor  axis.  As  one  traverses 
from  zero  elevation  to  90  degrees,  the  “closing"  of  the  fixed  polar  frame  at  the  pole  singu- 
larity causes  a fixed  error  in  dish  frame  to  subtend  a greater  bearing  angle.  This  behavior 
goes  as  the  inverse  cosine  of  elevation.  AH  we  have  done  then  is  to  substitute  for  o2  the 
value  o£d/cos2  E (subscript  D stands  for  dish):  the  inverse  cos2  E cancels  with  the  cos2  E 
that  appears  in  each  equation  with  o £ * and  the  subscript  D is  emitted  and  Equations  (9.17) 
are  obtained. 

While  maintaining  that  the  sensitivity  analysts  provides  the  ultimate  test  of  the  sub- 
optintality  of  a particular  filter  configuration,  it  was  also  felt  that  the  validity  of  the  error 
linearization  in  Equations  (9.14)  and  (9.16)  were  central  to  the  problem  and  that  further 
study  of  this  approximation  would  shed  light  on  the  sensitivity  results.  With  this  in  mind, 
two  simple  tests  related  to  the  linearization  approximation  were  conducted.  First,  it  was 
desired  to  directly  test  the  actual  assumption  of  error  linearity,  i.e..  Equation  (9.14).  This 
can  easily  he  done  by  calculating  the  true  nonlinear  Cartesian  error  as 


Vj-  = Xf  (xj>  + Vf)  - Xj.  (xp)  (9.18) 

where  the  function  xc  is  given  by  Equation  (9. 1 ).  The  error  due  to  nonlinearity  b then 
simply  the  difference 


SC  * Hcl  ~ vfc  (9.19) 

Actual  values  of  »&■  were  then  calculated  over  the  range  of  coordinates  of  interest  for  the 
C.FCS 
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which  exceeds  by  perhaps  an  order  of  magnitude  the  level  of  error  anticipated  with  the 
actual  tracking  sensors.  For  the  conditions  of  Equations  (9.20)  and  (9.21 ).  the  maximum 
difference  between  the  linear  error  and  the  nonlinear  error  never  exceeded  1 .5  yards,  i.e., 

!pi  < 1.5  yards  with  the  maximum  occurring  at  a range  of  35000  yards.  It  was  concluded 
that,  whereas  the  linearization  is  used  only  for  the  specification  of  the  measurement  error 
statistics  and  not  for  the  observations  themselves,  the  effect  of  this  approximation  on  sys- 
tem performance  would  be  negligible  since  (as  we  observed  in  a previous  section)  the  sensi- 
tivity to  utilizing  the  wrong  error  statistics  is  minimal.  It  must  also  be  kept  in  mind  that 
we  can  not  really  expect  to  be  able  to  estimate  the  sensor  polar  statistics  with  nearly  as 
much  accuracy  in  any  case. 

Another  interesting  test  was  conducted  to  evaluate  the  effect  of  the  nonlinearity  on 
the  transformed  Cartesian  distribution  functions.  In  this  test,  at  positions  over  the  range  of 
coordinates  of  (9.20).  one  thousand  points  each  of  polar  noise  were  generated  and  the 
corresponding  nonlinear  noise  values  vc  calculated  by  Equation  (9.18).  The  chi-squared 
test  was  utilized  to  evaluate  the  normality  of  the  errors.  Chi-squared  was  first  calculated  for 
the  original  polar  noise  Vp  (the  values  summed  for  the  three  coordinates)  and  then  for  the 
corresponding  sequence  of  vc.  The  results  were  at  first  very  surprising  since  for  every 
noise  sequence,  chi-squared  was  less  for  the  nonlinear  sequence  vc  than  for  the  original 
polar  sequence  Vp.  indicating  the  nonlinearly  transformed  Cartesian  noise  was  more  nearly 
Gaussian  than  the  original  polar  noise.  In  retrospect,  however,  the  reason  for  this  is  quite 
simple  in  that,  to  the  extent  the  linearity  of  Equation  (9.14)  is  valid,  the  Cartesian  erro.s 
are  simple  linear  combinations  of  the  polar  errors.  Of  course,  linear  transformations  of 
Gaussian  distributions  are  also  Gaussian  (see  for  example,  page  94  of  Bendat  and  Piersol 
(1971 )).  We  would  therefore  expect  by  the  Central  Limit  Theorem  (on  cit)  that  the  Cartesian 
errors,  being  approximately  linear  sums  of  the  polar  Gaussian  random  errors,  would  display  a 
tendency  to  be  Gaussian  also.  Indeed,  this  is  what  we  observe  from  the  chi*squared  tests. 

This  further  reinforces  the  author's  contention  that  the  linearization  provides  an  excellent 
approximation,  and  given  that  the  input  polar  observation  errors  are  Gaussian  the  corre- 
sponding Cartesian  errors  should  also  be  Gaussian.  (Note  that  it  is  possible  to  analytically 
describe  the  true  nonlinear  Cartesian  distribution  function  given  polar  Gaussian  noise  by 
using  the  Jacobian  of  the  transformation.  Unfortunately,  the  author  was  unable  to  simplify 
the  result  into  a form  that  is  recognizable  for  purposes  of  interpretation.) 


C.  UNCOUPLED  LINEARIZED  FILTER 

The  two  Cartesian  filters  presented  so  far  have  resulted  in  ninth-order  state  and  co- 
variance  systems.  Upon  examination  of  the  previous  filter  equations,  it  becomes  immediately 
obvious  that  the  only  terms  which  couple  the  three  directions  are  the  off-diagonal  term  of 
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the  Rcl  matrix  since  the  target  dynamics  were  assumed  uncoupled.  If  the  off-diagonal 
terms  of  Rcl  are  assumed  to  be  negligible  (an  assumption  to  be  tested  shortly),  then  the 
system  can  be  represented  as  three  third-order  systems  (or  “channels”  as  we  call  them). 

This  results  in  an  enormous  reduction  in  the  computational  burden  if  it  can  be  accomplished 
without  serious  degradation  of  performance.  In  fact,  the  computer  time  utilized  on  the 
CDC  6700  for  each  implementation  was  monitored  and  the  uncoupled  filter  found  to  be 
less  than  one  twentieth  the  other  two. 

Again,  it  is  interesting  to  directly  test  the  validity  of  this  assumption  so  that  we  might 
better  understand  the  results  of  the  sensitivity  analysis  to  be  presented.  A very  brief  con- 
sideration of  the  off-diagonal  measurement  covariances.  Equations  (9.17d-f),  is  very  in- 
structive. For  example,  each  cross  term  is  weighted  by  a factor  of  form 

a = Sin  X Cos  X 


= 0.5  Sin  2 X 


It  is  very  easy  to  stc  that  the  maximum  value  of  a is  0.5  and  its  average  magnitude  or  root- 
mean-square  is  only  0.  i 25.  Since  the  remaining  factors  are  no  greater  than  the  diagonal 
variances,  this  means  that  the  cross-correlation  coefficients  will  not  exceed  these  values. 
Further  examination  of  the  terms  in  parentheses  is  even  more  revealing.  Notice  that  these 
terms  always  algebraically  subtract,  thus  tending  to  further  reduce  the  correlation  coef- 
ficients. In  fact,  for  typical  values  of  polar  sensor  statistics,  the  cross  covariances  all  com- 
pletely vanish  at  a range 


rS  = °t  l°\ 


where  ox  is  the  angular  (bearing  or  elevation)  standard  derivation.  Depending  upon  the 
particular  sensor  (and  other  conditions),  r$  should  usually  be  between  1000  and  1 2000 
yards,  with  approximately  5000  yards  most  likely.  Clearly,  this  range  of  r$  coincides  with 
the  tracking  ranges  of  interest  for  Naval  gunfire  control.  For  all  these  reasons,  we  would 
expect  the  cross  correlation  of  channels  to  be  very  weak  and  the  suboptimal  uncoupled 
filter  to  be  close  to  optimal  in  performance. 


D.  SENSITIVITY;  TEST  CONDITIONS  AND  RESULTS 

In  order  to  evaluate  the  sensitivity  of  the  two  suboptimal  filter  configurations,  a set  of 
test  conditions  specifying  polar  track  sensor  statistics  and  targets  must  be  chosen.  Two 
hypothetical  track  sensors  were  parametrically  determined  in  Table  9. 1 A.  An  older 
"sloppier"  radar  called  Sensor  “A”  track*  to  only  3 milliradians  and  suffers  from  serial  cor- 
relation of  bandwidth  1 .88  Hertz  due  to  the  absence  of  off-axis  corrected  tracking.  A newer 
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Table  9. 1 A.  Track  Sensor  Statistics 


Table  9. 1 B.  Target  Characteristics 


MANEUVER  STATISTICS 


Maneuver  Frequency 
Nonmaneuvering  Target 
Maneuvering  Target  “B” 


wm  = 1/20  seconds'1 
oma  = 0-1  yards/sccond2 
Omb  = 5.0  yards/second2 


NOMINAL  TRAJECTORY 
Closing  Target 
Crossing  Target 


r = 17500  to  0 yards 
X = Constant  = 2500  yards 
Y = - 7500  + 300  • t 


VELOCITY 


Mach  Number  = 1.0 


TARGET  FLIGHT  TIME 


T = 50  seconds 


“tighter”  radar  called  Sensor  "E"  adds  one  milliradian  of  white  noise.  Both  sensors  track 
with  5 yards  of  white  noise  in  range.  The  targets  are  specified  by  choosing  random  maneu- 
vers with  given  statistics  superimposed  upon  a nominal  straight-line  target  with  given  speed 
and  orientation  with  respect  to  the  ship.  These  conditions  are  presented  in  Table  9. IB. 

Two  levels  of  were  chosen  to  represent  both  maneuvering  and  nonmaneuvering  targets. 
Two  nominal  transonic  target  paths  were  chosen:  oni  closing  directly  from  an  initial  range 
of  1 7500  yards:  and  one  that  crosses  with  a point  of  closest  approach  of  2500  yards.  Fifty 
seconds  of  target  track  data  were  assumed  to  be  supplied  for  estimation  and  prediction 
purposes. 

The  sensitivity  test  results  shown  in  Table  9.2  are  given  in  terms  of  the  prediction 
criterion  of  Equation  (9.8).  First,  consider  the  performance  of  the  coupled  linearized  filter 
relative  to  the  optimal  nonlinear  filter.  In  Section  9.B..  we  discussed  and  presented  some 
auxiliary  results  which  indicated  that  the  linearization  approximation  was  quite  accurate  for 
our  application.  If  this  linearization  were  perfect,  the  linearized  coupled  filter  would  be 
optimal.  We  find  that,  indeed,  using  our  criterion  of  optimality,  there  is  essentially  no  dif- 
ference in  the  performance  level  of  the  coupled  linearized  filter  relative  to  that  of  the 
optimal  nonlinear  filter.  We  also  note  that  the  performance  is  a strong  function  of  the 
sensor  statistics,  which  of  course  we  would  expect  since  the  accuracy  of  the  linearization  is 
a function  of  the  error  level  itself.  The  author  therefore  concludes  that  the  linearized  filter 
is  an  excellent  approximation  to  the  optimal  nonlinear  filter. 


Table  9.2.  Three-Dimensional  Sensitivity  Results 


Track  Sensor 

Maneuver  Level 

Target 

0 (Linearized/ 
Optimal) 

0 (Uncoupled/ 
Coupled 

A 

A (Low) 

Closing 

0.994 

0.956 

A 

A (Low) 

Crossing 

0.994 

0.961 

A 

B (High) 

Closing 

0.998 

0.978 

A 

B (High) 

Crossing 

0.9985 

0.993 

E 

A (Low) 

Closing 

0.99984 

0.978 

E 

A (Low) 

Crossing 

0.999996 

0.995 

E 

B (High) 

Closing 

0.9999998 

0.997 

E 

B (High) 

Crossing 

0 9999996 

0.9997 
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Now  let  us  consider  the  uncoupled  linearized  filter.  This  filter  is  most  important  since 
it  results  in  a significant  computational  advantage.  Since  the  linearized  coupled  filter  is  so 
close  to  optimal  for  our  problem,  we  compared  the  uncoupled  filter  to  it  to  measure 
optimality.  (Presumably,  one  might  estimate  the  criterion  for  this  filter  relative  to  the 
optimal  nonlinear  filter  (if  desired)  by  calculating  the  product  of  the  two  values  of  0 al- 
though the  author  docs  not  prove  this.)  We  find  that  the  drop  in  performance  is  more 
noticeable  for  this  case  but  remains  within  a few  percent  of  the  optimal  for  all  cases.  Again, 
for  the  improved  E track  sensor,  the  difference  is  essentially  negligible.  The  author  there- 
fore concludes  that,  due  to  the  significant  computational  advantage  enjoyed  by  the  un- 
coupled linearized  filter  and  to  the  relatively  small  loss  of  optimality  of  this  design,  that  it 
is  to  be  recommended  for  implementation.  It  is  also  felt  that  other  conditions,  such  as 
imperfect  target  motion  modeling  and  inaccurate  knowledge  of  sensor  error  statistics,  could 
result  in  more  significant  performance  reduction  from  the  optimal  performance  level.  This 
same  conclusion  has  been  reached  by  other  researchers  such  as  Cantrell  ( 1973)  and.  in  fact, 
independent  Cartesian  filtering  is  presently  employed  in  all  U.  S.  Navy  FCS  for  major- 
caliber  guns. 


E.  POLAR  AND  HYBRID  FORMS 

As  mentioned  previously,  since  the  uncoupled  Cartesian  filter  was  the  one  ultimately 
chosen  for  implementation,  the  details  of  work  on  the  other  forms  will  not  be  presented 
here.  A brief  discussion  of  the  polar  and  hybrid  forms  is  indicated,  however,  since  they  are 
both  interesting  to  compare  and  contrast  with  the  Cartesian  filters. 

The  spherical  polar  state  vector  is 

Xp3D  = lr  i t B B B E E EjT  |9.:2> 

The  kinetic  relationship  between  xp3D  and  can  be  obtained  by  doubly  differentiating 
the  nonlinear  transformation  Equations  (9.1 ) am*  (9.2).  The  dynamic  relationship  is  then 
determined  by  substitution  into  the  three-dimensional  Cartesian  target  motion  Equations 
i 9.5-7).  As  might  be  expected,  this  manipulation  is  quite  tedious.  The  resulting  polar 
dynamic  equations  are  highly  nonlinear  with  rather  strong  coupling  between  coordinates. 
The  state  extrapolation  can  be  effected  exactly  however. 

Unfortunately.  if  was  not  possible  for  the  author  or  his  associates  to  determine  the 
exact  nonlinear  extrapolation  for  (lie  error  covariance,  and  it  was  r.ecessaiy  ;n  linearize  and 
iterate  again  in  order  to  construct  the  optimal  nonlinear  polar  filter.  The  update  process  in 
polar  coordinates  is.  of  course,  linear  and  posed  no  problems.  Of  course,  as  far  as  perform- 
ance is  concerned,  the  optimal  nonlinear  filters  in  both  coordinate  systems  are  exactly 
equivalent.  If  one  were  determined  to  implement  an  optimal  nonlinear  filter,  then  the 
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choice  of  coordinate  system  would  be  based  upon  other  considerations  such  as  computational 
burden.  This  conclusion  is  essentially  identical  to  that  of  Monzingo  (1972)  who  decided 
that  both  approaches  could  be  made  to  yield  comparable  track  accuracies. 

The  basic  computational  problem  centers  upon  the  relative  degree  of  cross  coupling 
and  nonlinearity  (however  one  might  define  such  an  effect)  between  the  nonlinear  upda. 
step  for  the  Cartesian  filter  and  the  nonlinear  extrapolation  step  for  the  polar  filter.  This 
relative  effect  is  clearly  a function  of  the  particular  application  and  would  determine  the 
number  of  iterations  required  to  effect  convergence.  For  example,  for  ionp-iange  tracking 
of  reentry  vehicles.  Mehra  ( 1971 ) found  that  the  update  nonlinearity  was  more  severe  than 
the  dynamic  nonlinearity  for  the  assumed  conditions  of  that  estimation  problem.  This 
could  be  expected  since,  as  mentioned  previously,  the  observation  error  contours  are  very 
distorted  at  long  ranges  and  the  reentry  vehicle  dynamics  does  not  usually  involve  high 
angular  rates.  The  author  suspects  (but  did  not  demonstrate)  that  the  opposite  is  the  case 
for  our  application.  (It  should  be  noted  also  that  the  observations  with  the  phased-array 
radar  in  Mehra’s  paper  were  not  of  angles  but  of  direction  cosines  for  which  the  nonlinearity 
is  different.)  It  is  also  accepted  by  many,  including  the  author  and  Monzingo  (1972).  that 
the  independent  (uncoupled)  suboptimal  polar  filter  could  degrade  seriously  at  close  ranges 
due  to  the  neglect  of  the  "apparent  accelerations’'  mentioned  previously. 

The  concept  of  the  hybrid  filter  is  simply  an  attempt  to  minimize  the  effect  of  the 
noniinearities  and  to  use  the  best  (most  linear)  parts  from  both  the  Cartesian  and  polar 
filters.  The  state  vector  for  this  filter  is  mixed  (hence  the  term  "hybrid") 

in 30  = ir  B E x x y y z zlT  (9.23) 

The  only  nonlinearity  that  appears  for  this  filter  is  with  the  extrapolation  of  polar 
positions  since  the  other  state  variable  extrapolations  and  the  update  are  linear.  The  error 
covariance  matrix  is  still  full  (9X9'  and  must  be  linearized  for  the  extrapolation  step. 
Additional  work  with  this  filter  is  required  before  conclusions  can  be  made  but  it  remains 
an  interesting  concept.  Variations  on  this  filter  have  been  conceptually  explored  and  in- 
clude a rotating  Cartesian  filter  with  origin  maintained  at  the  current  estimated  target  posi- 
tion and  orientation  such  that  one  axis  is  aligned  with  the  ray  and  the  orthogonal  axis 
aligned  tangent  to  the  sphere  in  the  directions  of  bearing  and  elevation.  Clearly,  additional 
investigations  are  needed  if  one  is  to  make  definite  conclusions  as  to  the  advantages  (if  eny) 
of  these  ideas. 


X.  SUMMARY  AND  RECOMMENDATIONS 

This-  section  is  intended  to  give  the  reader  a brief  overview  of  the  studies,  results,  con- 
clusions. and  areas  for  future  work  contained  in  this  report.  Section  I presents  the  reader  a 
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brief  introduction  to  the  fire  control  problem  and  the  basic  considerations  attendant  to 
"front-end”  (filter/predictor)  design.  Section  II  introduces  the  concepts,  notation  and 
equations  of  the  Kalman  filter  with  particular  emphasis  on  reasons  for  its  choice  in  this 
application  and  information  required  to  construct  a filter  of  this  type. 

The  next  six  Sections  (III-VI1I)  deal  exclusively  with  modeling  and  parametric  be- 
havior for  one-dimensiona.'  filtering  and  prediction.  The  derivative  r ‘lynomial  target  model 
was  first  explored  (Section  Hi)  since  it  is  undoubtedly  the  most  common  and  familial  and 
parametric  behavior  of  a rather  fundamental  nature  was  unavailable  for  this  simple  model. 
Convergence  studies  indicated  that  constant  jerque  and  higher-order  niters  could  not  be 
considered  for  this  application-at  least  as  the  principal  fiiter-since  they  take  an  unaccept- 
able amo**nt  of  time  to  adequately  settle.  Quantitative  relationships  of  filter  convergence 
with  polynomial  order  and  observation  data  rate  were  determined.  The  problem  of  dealing 
with  a maneuvering  target  with  unknown  strategy  led  us  to  consider  the  random  acceleration 
target  model  in  Section  IV.  The  author  found  this  model  particularly  appealing  since  it 
statistically  recognizes  and  acknowledges  this  lack  of  information  and  enables  one  to  directly 
rclat-.*  the  targe:  maneuver  parameters  (acceleration  level  and  fiequency)  to  the  Kalman 
filter  bandwidth  and  performance  matrices.  It  was  demonstrated  that  this  filter  converges 
faster  than  the  const  an;  acceleration  filter  and.  in  fact,  tends  to  display  behavior  somewhere 
between  the  constant-acceleration  and  constant-velocity  filters.  This  model  is  particularly 
valuable  in  that  it  enables  the  user  to  determine  fundamental  limitations  concerning  the 
ability  to  filter  and  preuict  a maneuvering  target  path.  Additional  work  on  target  motion 
analysis,  especially  from  a *hree-diniensional  point  of  view,  is  recommended. 

Concepts  of  divergence  prevention  and  adaptive  Kalman  filtering  arc  considered  in 
Section  V.  A discussion  of  divergence  and  the  bandwidth  tradeoff  is  followed  by  the  devel- 
opment and  analysis  of  the  urw  of  residua;  analysis  as  a maneuver  detection  tool.  Having 
a divergence  detection  technique  in  hand,  two  concepts  of  adaptation  are  developed.  The 
advantage  of  the  dual-bandwidth  over  the  single-variable-bandwidth  adaptive  filter  are 
presented  along  with  the  particular  adaptive  scheme  recommended  for  implementation.  The 
dual-bandwidth  concept  is  particularly  attractive  since  both  narrow  and  wide  bandwidth 
filters  arc  always  in  operation  and  the  adaptation  algorithm  merely  selects  the  output  state 
vector  according  to  the  monitored  performance  of  eaci:.  Additional  work  relating  to 
advanced  adap’ation  techniques  such  as  multi-bandwidth  parallel  or  cascade  filter  banks 
with  hypothesis  testing  or  residual  analysis  and  in- need  *>•  tsjrivity  adaptive  designs  are 
strongly  recommended. 

The  problem  of . erially  correlated  observation  cr  ors  is  treated  in  Section  VI.  A cor- 
relation model  ami  modified  Kalman  filter  and  sensitivity  algorithm  arc  developed,  and 
results  pertaining  to  observation  error  .statistics  are  presented. 


It  would  probably  not  have  seen  possible  to  actually  implement  these  Kalman  filter 
designs  on  a real-time,  fixed-point  ^uni-computer  if  it  were  not  for  the  prefilter  and  square- 
root  covariance  techniques  developed  £•>  Section  Vil.  The  data  compression  work  in  partic- 
ular is  critical,  and  more  advanced  state-of-the-art  techniques  should  be  explored  in  the 
future  since  the  payoff  in  reduced  computation  is  so  great  with  a successful  design. 

The  predictor,  of  course,  is  the  ultimate  product  o <-  cm  front-end  (target  sensor 
plus  the  filter)  so  a special  section  (VIII)  was  dedu^sc  *«  of  predictor  design. 

After  demonstrating  that  the  minimum  variance  preu.t  r oy  the  transition  matrix, 
the  sensitivity  of  the  predictor  to  the  choice  of  accelerator,  -me  constant  and  prediction 
accuracy  for  the  limiting  case  of  perfect  filtering  was  discus  *d. 

Finally,  in  Section  IX.  the  results  for  the  one-dimer*  ' 'her  work  were  tied  to- 
gether into  the  three-dimensional  filter  design.  After  devei,  ?>ng  the  optimal  nonlinear  filter 
as  a standard,  it  was  shown  that  the  linearized  observation  error  covariance  approximation 
is  quite  accurate  for  our  purposes,  and  that,  by  paying  only  a small  cost  in  performance 
degradation  by  neglecting  the  cross-correlation  terms,  the  three-dimensional  filter  can  be 
decoupled  with  a dramatic  reduction  in  required  computational  load.  It  is  recommended 
that  additional  work  on  hybrid  and  polar  filters  be  performed  along  with  the  target  motion 
analysis  problem  with  the  objective  of  further  performance  improvement  foremost  in  mind. 

The  final  version  of  the  FORTRAN  mode!  incorporating  all  the  features  discussed  in 
this  report  is  found  in  the  Appendix. 

While  it  is  felt  that  many  improvements,  even  of  a fundamental  nature,  can  be  r.iade  to 
this  filter  design,  this  model  represents  the  best  that  the  author,  with  the-  aid  and  advice  of 
many  associates,  could  achieve  to  date.  It  b believed-  or  at  leas!  fervently  hoped-that  this 
model  represents  an  advance  over  existing  filter/pred  jetors  currently  employed  for  weapons 
control  purposes. 


106 


SUBROUTINE  0fclVti<(Ly8XVf6YVy UZV) 

GUM  MON  /FILTER/  0T»1EMf£LRH<MtXLPH4ifi«£TiA#i!'T4Sf  GAKFABf 

2 aAUyQt>ll.UAl2ylM12yaA13»U3i3yCA22ya322yOA2JyQ;>23yCA33yeoJ3yGly 

3 GZfCZtiTO^.ex  »S r0»t£2»ST0Ri3f MtHZ»f'£«A3f  MJk33» SIG3NA»SIG3N8»SI62*#» 

4 HC23»RH023,  Fhf5Ny  T AUPyF  1,  F2  y0ATOT£.  C.  PO  lNTSySIGHAMA  ,S;Gh*M8» 

5 TAiHA,TAU-f0fRHO12fSlGM>2»SlGR2»3XG62  ,SIG -2 
REAL  f£M?,*:MAJ,a'N33 

COMMON  /STATE/  X4l*XA2yX£3,  Y£lyVA2yY*3yZ-il,ZA2*Zl3yUUMrtY£  (3>, 

2 *3i,Xd2, X93,  YBl,To2vYcG,Zol»Za2f  Zt>3rOUHMYi(3)  t 

3 SAXllfSAXtZf  G£X13y3AX22y5£  X27yS4X33y5AYily  :UY12yS£Y13y.»AY22y 

4 Str23,iar33,SiZlikS*Zl2,S^Z13,S4Z22,S*.Z?J,54Z33,ibxll,i«Xl2t 

5 33X13»SdX22»S&X23>»8X33*StfYU»S3Y12v SSY 1 3($tfY22,SBY2Jf  S3Y33y 
c S0ZU»S0Zi2»3bZi3yS8Z22t3tfZ23vS3Z3i 

P.IAL  UU4Xd~»» .^UAYOAEyNU»Zi*ARyNU3X3A3y.‘«U3V:i-RyNU3ZB£R 
G AVERAGE  MEASUREMENT j 

IF  (L  .EQ.  1)  GU  TO  12 
ICOUnT  - ICOUNTH 
zx  = zx  ♦ ?xy;f;i£m 
ZY  - Zf  * =Y tf/FMEN 
ZZ  5 ZZ  ♦ BZy/FmEM 
IF  t ICOUNT  .J.  t)  RETURN 
*3  ^OfiTiNUE 

C SPECIFY  SENSOR  STATISTICS 

CAll  CARTFUL  <ZX,ZY»ZZ»*,3»l:> 
m2  * ***2 

5IN32  = <SiM<e»i**2 
CUS02  * 1. 'SINGE 
SIN£2  = IS1M£I>**2 
C0SE2  * 1.-SINE2 
C PR-FlLTik  RtOjCTlCNS 
SIGRZ  * SIGPZ/Fnck 
SI632  = SXG82/Fk£k 
iIG€2  * SICEZ/Fr.E* 

C CARTESIAN  STATISTICS 

SIGX2  * SIN82*CQSE2*SIGR2*R2*CCSS2*S1C62 

♦ R2*SINB2*SIN'2*SIGE2 
5IUY2  * COiS2*CUS£2^SlGR2»R2*SlN32#SnB2 

2 *R2*C0SB2*SlN£2«SIG€2 

SIGZ2  » SIN;  2*SiGR2»R2*G0$i.2*SI&;2 
C CAwL  INOIPENOsNT  GAF.TESIAA  F.LUftS 
IF  CL. Ed,  1)  GG  TO  1 

GALL  FXLTcR<X«l»XA2fXA3tS4XUf$£Xl2,S*X13»SAX22,S«X23,S*.X33» 

3 ZXfSIGXJ,  N0AXfc,.RyNU8XB*iR) 

1 X61»X4i2vX63ySGXlltSuXi2ySUXi3tS3X22fS3X23ySax33y 

2 ?0itYB2yY43fSBYll,SBYl2y5dV13»SSY22,33V23'S4V33y 

CALL  Fit TiP.LYAly VA2yYA3»SAYllySAYl2t> AY  i 3y>AY22yd»Y23yS~Y33» 

2 ZYyjIGY2  yNt)A Yt»ARyNUBYBARJ 

2 ZBlf ZB2tZ93ySBZllyS6Zl2fSBZ13y$3Z22y SBZ23ySBZ33| 

CALL  FItTtK«ZAlyZx2yZi3ySAZllySAZi2i,3iZiJy3:Z2£y3AZ23ySAZ33y 

3 ZZySXGZ2yNUAZt>»*yNU3ZB4R> 
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iCOUtT  - 

zx  = : 

2T  s G 

22  = a 


: SQr.T (S  Iu*2) 

ISIGY  = 

*QftT<SIt 

- ST  3+.i  3* ilGx 

HcXIi 

* S.xn 

- ST Crt  3*  GiG V 

SS-.Yil 

* SAYii 

= STy*c3»»lGZ 

3^b2ii 

* Stzii 

* *T’3CL2* jiGX 

SSjX  < 2 

- S-X12 

= ST0Ki2*aiG» 

3S-JY12 

- s;»i2 

- STykt2*jI62 

CSeZir 

* S4212 

- ^HC-i3*SiGx 

ISBxid 

<4X13 

i;  run* 

i Cu.'«T.HU= 
c STATE 

Xsi  = dXY  >Xcs:  S X A 1 

= 5t¥  jrei  = fti 

1 = SZx  iZtii  = 2ii 

C £R30<%  CjXAMiOCc  SCUAkc 
Si  OX  = 

Saxu 

s«yi; 
sazu 
Sixtr 
*4*12 
SA  2 12 
s*xis 
satu 
SAZU 
$4X22 
SAY  22 
S4ZZ2 
SAX23 
S4Y23 
SA223 

S»X33  - »*G3MA 
SA»23  = SiiJSA 
S-223  s SIG3N4 
NUAxiiR  s 
NUArjhfc  x 
MUA23Ai>  * 

icou:«r  s ; 

2X  s c. 

ZY  = C. 

22  = e. 


3X4  2 
3V-2 

*22  2 
OCT 


* o. 


3X3* 

3X52 

3232 


**»»OX3*$IG»  3S6X1S  x SiU3 
^rtui3*S-G2  3Sc2i3  s S4213 
5Ty*£l*5IG2M  353X22  - 54X22 

S2K22  3 53  X 22  * 54X22 

S4*22  I5SZ22  - S4*22 
^Hj23*SiG2N  3S5X23  * $£* 23 
32X23  35dY23  = 34X23 

SxX23  SS3223  = $4X23 


356X33 
iSc Y33 
356,233 
ifiUjxbAA 
tOUSYSAK 
3.XUS23AR 


SlG3-'<8 

SIG3N6 

SlG3n& 

** 


12  CONTINUE 

C CALCULATED  CONSTANTS  - REJUifcCU  BY  FIlTE* 
UT  * FtoXTcm N)/OAT«ATE 
gahnaa  = :xP(~oT/rAuHA) 

BtT-.-.  = T-iUNA*  <1  .-GAHHAi) 

G A lined  s r XP(  •DT/TAUtt3) 

BiTAB  = TAUKbMi.-GANNAB) 

ALPHAA  = UUflA*‘2»  cGA*tUi*DT/TAl*it-l.» 
AlPmAB  = T4UtB**2*(GAHHC3»0i/TAU!l8>l.3 
QA33  s 2«*yT*  SiGH4HA**2/TAUf1A 
<3833  * 2. •0T*SlGHAH8*»2/TAUHi 


s - 3X43  = 3X63  * 

* 3*  3Y„ 3 s *,  Jrb3  * 

1 »•  3243  * w.  $263  * C. 

3aiS2  = SaRT(SIG22l 


I OS 


♦ 

It 


4 


- 01*^33/2. 

Qd?3  = 0T*«J83  3/2» 

Q*£2  - 2,  -DT* QA23/3 • 

3gz 2 ~ Z.*3T*Q623/3. 

* UA22/2. 

<*•-'  3 = CS22/2. 

3-12  = i.  •0T*Q<.l3/*». 

Q3i2  = 3.  •0T*Qf»l3/G. 
a*U  * 2.*07*0il2/&. 
asu  = 2.*or*Qbi2/s . 

&2  = i./PClNTS 
= 1.-G2 
£2  - G2*C**2 

iTO^il  = SORT  Cl  -*h023**2) 

^rna-f  = <*HG12-*rtO13*RH023>/STO*ti 

K£*  " CX‘*S1^ ME2-2-*H0l3?.2i 

P«£h2  * Frl£ i»« *2 

SOI  = C. 

SOMA  * C. 

SUM6  = C. 

03  230  I s It 
rl  = I 

SOM  S SUMffi 

SUK5  * S0'H*GiK.»AA**CFl/FHE<o 

*-c 

»OM  = SU«/Frul»C2 

SUM-  s £UN./(Fr.ti«,6iH!Ui| 

aOiid  * SO*d/  f FRc4*GAIImAd) 

M-M2  * OTMl.-Sorf) 

* ;'',JlHf**2‘Csy',"-i**-rAOM£*C;i#2 

r*  - CFMcKfI#  )/(2>*FH£|i) 

F2  = tftlik-l.  )/(2#  •f«£„| 

SI&3JC6  s SXGMAK- 
SIG3K8  = SIGMAKi 
Z*  * ex* 

z»  * er* 

ZZ  = 82* 

GO  TO  13 
£*»0 
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uo  ouoooii  o nuooo 


C 
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SUBROUTINE  FicT. t(X Al, X»2,  X .3,SA 11 ,3i 12,S-12 ,3£Zt , SA£ 3,CA33 , 

2 Xol,XSZjXL2,Sbll,Sdi2,S 313,5822,552?  ,So33, 

3 ZySIGHA *2,NU At- R,hUrJiAP> 

LAST  UPOATE  - 2L  JOLT  197*. 


QNi-OIHINSIONAl*  ADAPTIVE  GUAL-3AN0WIDTH,  KA.HAH  TRACKING  FILTER 
Filter  a - narr.w  e-ndwioth 
FILTER  t i *IuE  dANOnIOTM 


VARIAoL-5 
XA/X6  * 

1 = 

2 = 

3 =■ 

SA/Sd  • 

l - 

3IGHLN2  s 

CONSTANTS 
OT 


STATE  VECTORS 
POSITION 
¥£ LCCiTt 
ACCELERATION 

Er.ftJ<E  COVARIANCE  SAOirt  ROOT  liTnlCii 
NcAiuRENtNT  (AVERAGE  OF  KE*  335ER*~T IONS) 
VARIANCE  OF  WHITE  «£A SURESENT  £*<CR 

TlH£  INCR-ML NT  dlTNE-N  FILTER  CTC_ES 


CONNON  /FILTER/  DT * rt€N»  ALPHfcAf  AtPhAg,  BET  AA,6cTAa,G  A*NA*  y&J-HHAd, 
c Q£ll,Q511,Gtlt,Q9l2,0<<13,G3l3,34  22, 2322, 0423,0323, Qi 53, 0853, Gl, 

3 G2,C2»!'T3N£*fSTC4U2,STcR£3,M£tt2»r£4A3v'i£NB3,5iG3NA,SlG3Ne,SIG2N, 
« H023,*irt023,Fh.M,T»U*',Fl,F2,0ATRfcT:.:.P0:NTS,Sl&NArti,SiGrA«0, 

£ TA(MA,TAL?1&,RNG12,SlSn.2fSl£r'2,SiG42,SlG£2 
RE  At.  K>1,<m2,  KA3,R»  l,<32,<di,NUA  ,NU3 
RiAt  NFN2,''-NA3,*irN32  ,NUAa«R,NUd3~R 
STATE  EX-RAPOwiTI<Mt 


X*i  = XAi »UT*XA2»AlPHA£*XA3 
Xfll  * Xo I *0 1 *Xc; 2 VAt. PH 43*  XB3 
X-2  s XA2*3.T£A*XA3 
X82  * XB2*8ET A8* X33 
XA3  = GAH:3fcA*XA3 
Xb3  = GaNNA  3*  Xbi 

SOU-KE  SOOT  COVopvIAkCE  EXTRAPOLATION 

TIGH'aO  = saRTCSAil**2*SAl2**2*Sfti3**2*$lGJ!£a2) 

ifc*  .J60  * SQ*TCS3il**2*S3t2**2*S8l3**2*i:CfU*2) 

SA  12  = SA«2*0T*$A  22 

SM  i Sai2*OT*$322 

Z*>  * S*l3*0T*  SA23»AkPnAA*  SA33 

io.  . * Sdl **U'»S323*ALPHA6*SB33 


3£-i  s 5ft2J*&ETiA*S*35 

SCtl  - S023*8£T  A3* S 83 3 

CA-TO  • GAH-aA*SA33 

S3  2 A GAHHA8*Sd33 

*833  S 3833* SORT (1*  *Q333/S333**2) 

Ju23  * <SdJ3/*t*33>*Sd23*Qe23/w633 

*922  = S822*S0R7U.  *lf$823**B23)*  ($323-*i2J)*0l*22)/5B22-*2> 
<*81 3 * (S  E33/M833)  * 5813*091  3/HB33 

N«12  - (S822/ *£22) • S812*  <S6 13*S623»*B1 3**823*  Q812) /MB22 
*811  * *811*SQRT<i.*«IS8i2*w8l2)*|S812-W3l2)*CSSl3**8l3># 

2 ( S3l3*M8l3) *0611) /SO  11**2) 

CORRECT  nEASURENENT  FOR  TARGET  nOTlCN  C,HEN  P*£FI LTERING) 


no 


i 
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c CAwcucrt  a Pkiohi  f.'sxoual 

.<U~  = Z--x  :i 
*U3  * Zd-Xtji 
c «».f»wuvis  ofiEicroK 

s*»  = PAil*SI6HAM2  * 

* Pen 

5»G*UA  = SQKT (&„ ) 

S.G?<U6  i SdP.T  iQit> 

SwiT-i  I 

r *i£»  0*>  fiO  TO  1 

C Mhcjitn.  u£C  L~rtATiQ<<i 

= xai 
*»2  s <d2 
***i  = XB3 
*0«  = NUd 

U‘tf6fc  * <S£6N0o/SX0MU^»*MUd84K 

*?f*  * f:JJ‘SGKTli.  +OA33/.S4 3J»*2» 

* tSi3?7^35,1;sJ{|^|jJ^fJJ<S*zJ-»-2J>*««22»/Sil22**2> 

•^il  * S*ii*S»TU,*<cIli|*^5jf3ZffZ'3*l'AZ3*Qa12»/»»*22 

2 (SAU-kAllUQSui/SAip^)  ‘ sl12*WlZ>MS»13»MiJj. 
P-U  * P<tXl*(,0lx*U3ii 
* CA-tUlifQfrU 
GO  TO  «. 
c HO  <UN£U«E>( 

1 COHTlMOi 

J*f3  * X??tSaftT  <l* 

£ff  - shJ2J^i,,*s*M#®t2,/Mw 

**z‘  * $<22 *S(*»FI  !•  *f  # 

;j{j  * !!:«,,wj3»*seiSQ*i5/5?33  2J,*°4tt,/S4?2**?» 

Xin  * slijL*s  OPT fi Tf ?clr If i iJIf ff 

2 (SA13~mAX4)*QA££)/SA  ••*  #,»  ^ 1 2*fc'UZI»  (SXtltHtx})* 

8 COMTlMUt  44  ** 

0 GAINS 

!?}*  * *f»**«3Z2M«8i3*Ntt2J 

P#»1J  * NAX3*NA33 

P»13  * *613**833 

*»i  * P*n/GA 

<81  * P«tll/G6 

<A2  * Pilt/Gi 


m 


t 


K42  * PB12/GB 
Km3  ' PA13/GA 
<63  * PdlJ/G8 
C UPDATE  STATE 

FACTOR  = i./  IFi*f2*SIGNUiO/Sl6NUA) 

FACTuRB  = 1./ CFi*F2*SISNOdO/Sl6NUei 
KuA  * F«.C70R£**U* 

NUB  = FACTORB’NyB 
X*  1 = X*1+<U1*M)a 
XBl  = X£l*K8l«NU3 
X<*2  = XA2»<*2*KUA 
XB2  = Xe2*«82*Nua 
XA3  = XA3»fc£3*hU» 

Xt»3  = X83t*B3*)»Ud 
C UPDATE  GOtfAAUWCE  SQUARE  ROOT 

SA33  = NA33*iQRT (1* -6A*K«3*XA3/MA33* * 2) 

S333  * nd33*SQRT fl * “B8*K83*f(B3/MS33**  21 
SA23  * <NA33/S*33)»MA23>6**KA2*K«3/SA33 
S623  * (N333/S833I*  frB23~GB* KB2*K33/S333 

SA 22  = NA22*sDRT  <1,  ♦<(NA23*SA23>*INA23-SA23I-«A*«A2*ICA2I/MA22#*2» 
SB22  * NB22*S GF.T (1  • ♦ ( (NB23*-SB23>  *<N823-SB23l-CB#K82*Kl»2)/Nfl22##2) 
SAA3  * (NA33/S*33>*MAl3'tA*<Al*KA3/S03 
S813  * (N333/SB33I *M6l3-GB*Kil*KB3/5i33 

SA12  * iMA22s$A22>»MA12*<NAU»H&Z3-$A13»SA23~6A*KAl*KA2>/SA22 
SB 12  * INB22/SB£2)*MB12A(NB13*HB23>3Bl3»SB23>GB#KBi*KB2)/St22 
SA11  * WAlt»SQftTCl.»t<MA12»SA12>*{NA12~SAl2>»tlUl3»$Al3>» 

2 CNA1 3-SA13)-6A*iCti#KAl)/tiAil#*2) 

$811  * xiU*SQFT(l. W ()«12»SI12I* (N6l2*S|t2»M«313»Sbl3l • 

2 twai3-sei3i-ce**8i*<eii/NBii»*2i 
RETURN 
END 


Su^.aUTlN.  FAiOISf  (XfTP» 

CQNNDN  /FILTER/  OT,  .-1£ N,iLPH- A,  U»NAB,  BETAA*  AETa it SAMH&& , GihrtM, 

2 *iUfQili.w<.i2*33l2»  JA:3,«-»l3*CA2c,iJ22,C,,2J»QJ23,UA33*Q833,£i* 

3 »»2,32, S»GR-*  l,Sr3Rt2,ST0RE3»«EN2,*1£NA3,-ENd3»SI53NAfSlC3N6,SlG2N, 

«.  H023»xK023*  Ffi.  N,TAuP#Fi,F2,OlT<ATi.  C.  HJiNTSfSI&HAHA.SIGNA «B, 

b T-UnA,  TiUlU,rthOi2,:>lSHt2,SiW»S:..a2  ,SIGt2. 

DIMENSION  X (12) 

AlPH„P  S UuP**2*  (fc Xp  1-TP/TiuP)  *TP/TtUP-l .) 

X(l«»  = Xf  II  ♦TP**(2)**LPH*P*X(3> 

X(il)  * X(H*TF*X(*| ♦AlPRAP*X($» 

X(12)  s X(7|»TP*X(e)»ALPHlP*X(S) 

re  turn 
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listings  of  Computer  Programs 
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KALMAN  --  T?'i>  subroutine  perform?  ih*  »n»v?x  algebra  involved  in  one  cycle  of  the  general 
Kalman  niter,  The  general  Kahn?*)  iV'-r  assumes  s compicfe'y  linear  system  with  un- 
co rrelated  measurement  and  process  noise.  The  equations  sre  suwmsrizrd  in  Table  2.  i 
uiut  seproduceo  here  in  the  order  in  which  they  are  use/}  in  the  subroutine.  Noiice  that 
we  have  included  a delerminist'c  forcing  function  in  the  slate  equation.  This  inclusion 
nas  no  effect  on  the  renaming  equations. 


Xfk/k 

M - <j{k.  k - l X(k  - 1/k  ~ 

!)  + Ok  - )) 

(A,  1) 

P(k/k  - 

5)  » Ofk.k*  Ijl'tk  - 1/k  - 

\ rtr<k,  k - n 

+ Q(fc-  !} 

(A. 2) 

K*Si}  — 

P(k/k  - l)HTfk){H<k)P  fc/k 

- h»T(k>  + 

R(k>!  5 

(A.  3) 

Xfk/k) 

* Xfk/k  - l>  + K;/:>{Z(k> 

' H(k>X{k/k  - 

•Di 

(A.4) 

Pik/k; 

- H - ,Ltk>H(k>| t*<k/k - 1) 

(A.S) 

WW.  *tNc  >IH> 

u»R«?r  i>*  $f 1173 

r,FMFP»t  K»t*»in  rytree  cva.*-  rm*  u*»r»p  systfh  :i:tm  vimitc  kf '*'hfht 
bn*  ponces  intsK.  given  input  *r  tthe<k-i»,  this  n3jti<f  returns 

ftllTPtlf  *T  Tr*M>C». 

c * - sriTf  v-ern*  - tuouriK-ti 

r out  - TPimtnoN  mrs-.x • - input 

C nrrrpUTSTSTfC  EO*CTNr,  VECTO*  - INPUT 

C 9 - r»*0®  rf?V«»l*Nnr  X*?*T«  - I NP’JT  f < — | > - 9UT«»0TU» 

r.  1 - Nnt«;f  N*TPIX  - TN»UT 

CP-  Hc •SUPf  ^NT  *IOTSF  N6TPTX  - INPUT 

Ci-  nB^py-TTnn  h*t«ix  - input 
r K - r.BTN  «*TPTV  - C*lCUt*TFf) 
r 7 - »rasilPF,/rKT  yFCTOP  - INPUT 
C • .I.C.'’  - unvxTNC,  STOP*Gf 

r.  «m*,**nTN  - unfnmins  nF  vaPtBPirs  in  hailing  poutine 

r N,«*  - oorpftTfvc,  TT**rffSinNS 

C NITr  - f,  Pf  OT^FNST (INEQ  1PM  F03  THF  MATPIX  INVERSION  TO  <0S< 

PF*t  <C 

*?**«■«*  T nv  x:e'STM*.PMI{*jnrH,NTI‘H,r<wi«l<P(NTP1,NOIMI,U{NOINtNOtH 
Z ,Pf*nTN.H')IM).<<»NOTH,N')IH>  ,7INriN»,Hff*DrHfM0IH».4<N3IMI, 

T -MNOTISNnTNI  fCfMni*<tHriIP|,f)l'OTH,*1f)TN| 
t rxTP»pnt*TnN 

r s * ®Mr{K,tc-ti*»tx-i/x-i>*E«-i» 

C p r P>1  )*P|X-J/K-tl*TP*NSP3SrCP.-in*QtK-i) 

n»  t I * t,N 
»<T»  * Fff| 

’’I  t J * t,N 
• in  = - *n»oMT»r,ji»xt ji 
«<t,ji  - -H,J» 

mi  XV  r I,N 

on  t l * |,N 

nri,j>  = pii.j»*phht,«i*phi(j,{.i*p<xk,l> 

I HUNT i *mr 


A-l 


C <<*/«(-<»  * 4 

r PtK/K-t*  i o 

no  ? t = 

*(TI  * M!t 

n n ■»  j - t,*i 

■>  cnurynuc 

r 'jorjute 

r f.  s H(Ki  *Pf*C/<*-lt  • TP  f MfKI  I ♦ P 1 1(1 

my  T = t,« 

no  y J r 1,* 

rn.jt  t o<T,J» 

**')  i x<  = l,s 

•*0  ? L - i,  N 

r»T,jy  r 
t roM*  J»kic 

on*  p{e/r->«  •^*4*JSoo5f{H«tcn 
on  i r = t,* 

"Oft  J s t,f 

0«f,J»  = 0 

•0  ' V KK  = I , H 

* on, j» 

< ,rt*»Tt»uir 

IM*  .Cv.  11  OStt  *4TotX<tO»".*,3,r..“  ofT) 
if i*  cfi.ii  = i./cti.n 
o kik»  -•  o*,- 

oo  s r * t»*i 

no  ^ j = i.«* 

- o. 

•'0  & fc<  S 

* KtT,  Ji  ♦01T»«)«CS'0<*JI 
«;  rmif ‘Nr 

f ''(T,  1)  * m :-*■*»<»  »y  »*/<-!> 

».  i - ?»" 
r:r,i>  r mm 

•'O  * J : 1,1 

C * T . I > r .;:•*!  J) 

f.  r "ftt  t ».t!r 

*■  <<*/<!  = y fH/'C-i  UKIO  » 

"•or  r y i.n 

"o  y J S {,« 

ytTt  r )T(M»viTrJt*ou«l> 

■» 

F 0 : < «1  *M<  » *o»  K/ic-t  > 

"O  ' T I 1,1 
"o  a J * 1 , V 

n(y,j)  = a. 

"O  * n r |,y 
oo  a l - t , N 

PSI.J1  = Pn,J>»Ktt,<KI«vHKK,t»*0{l,J> 

* f TV1  T *,»)y 

r ofv/irj  , pfv/r-tj-o 
ci  1 t - t,<< 

oo  1 J = ;,*J 

° ! T * J I * -"a, 

Mtjc 

or 
r *jn 


RANDOM  This  subroutine  generates  a Gaussian  random  number  with  zero  mean  and  unit 
standard  deviation.  It  uses  uniformly  distributed  random  numbers,  generated  by  the 
CDC  system  routine  RANF,  to  compute  a normally  distributed  number. 

CORNUM  - This  subroutine  generates  the  ith  value  of  a sequence  of  exponentially  autocor- 
related  Gaussian  random  numbers  with  zero  mean.  The  output  of  this  algorithm, 
modeled  as  a first-order  Markov  process,  is  equivalent  to  a low-pass  filter  with  cutoff 
frequency  = l/(2nr)  Hertz  driven  by  white  Gaussian  noise.  The  sampling  frequency  is 
automatically  adjusted  to  satisfy  the  Shannon  sampling  theorem.  Additional  details 
and  documentation  of  a similar  algorithm  can  be  found  in  Ualas  <T967). 


subroutine  ^t*JDo*ttRAMNy-,:5 

c this  season.-*  c-:f..«-Tds  * randur  number  <rannu?o  from  a giussian 

C DISTRIBUTION  with  zero  MEAN  mo  UNIT  STANDARD  UEVJATIOT 
C sort  - IN  O-tDEO  TO  yOTAIN  A DIFFERENT  SIQINCE  CF  RANDOM  NUHJERS  FOR 
C tAC.:  CONFUTE.-.  RON,  R*hF  CAN  uE  RAN30«LV  INITIALIZED  WITH  ThE  SVSTEn 

c cluc*  av  s.rrr.x,  : * i the  first  Tint  this  routine  is  called 

C IF  lilt  SA**t  St  OuEl.cC  is  DESIRE 0,  St  T I s C THE  FIRST  Tit* 

C OTHERWISE  'ARE  1 .Nt.  1 
QUA  TwOPI/fc.  Zfc  J1553/ 

IF  «I  . GT.  1)  GO  TO  2 
IF  (I  .ED.  Cl  GO  TO  1 
SIGN  s Z7777 777/777/7777777* 

CtOCR  = T X* - E IDVhNT > .AND.  SIGN 
CALL  RAuScT tCLCCK) 

GO  TO  2 
i CONTINUE 

U..LL  F*NSrT(C.J 

: conjinu 

XI  = RANF(UUKHV) 

x?  — RANFtntnrj 

RANNUM  ; S*IPT  C -2«*A lOG( * II  I *CG5t  TM0PI *X  21 

RETURN 

ENO 
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SUBROUTINE  C03MM(XN0ISE  tSIGHA  t TAU»OElTAT  t I»N  »RhC»SIGHAH) 

c this  subroutine  generates  the  i-th  value  of  a sequence  of  auto- 

C CORREcATEO  RAMOGH  NUMBERS  WITH  ZiRO  HEAR 

C EQUIVALENT  TO  ; LOh-PASS  FILTER  NITH  CUTOFF  FREQUENCY  * 1/ (2*PI»TAU> 

C HERTZ  DRIVEN  i»  MhITE  G.'USSIAN  NOISE 

C XN01SE  - OUTPUT  HE*  RANDOM  NWBER  - INPUT  OLD  VALUE 

C SIGH*  - jTANOA-'D  DEVIATION  OF  OUTPUT  SEQUENCE 

C TAU  - EXPONENTIAL  AUTOCORRELATION  TINE  CONSTANT 

C OELTAT  - CONSTANT  SAMPLING  INTERVAL 

Cl-  INDEX  OF  SEQUENCE  < .GE.  1) 

C THE  FOLLOKING  QUANTITIES  AP.f  COMPUTED  INITIALLY  (1*1)  ANO  MUST  SC 
C STORED  IN  THE  CALLING  ROUTINE  *0*  USE  HMEN  I .GT . 1 
C N - FACTOR  Of  INCREASE  OF  EFFECTIVE  SAMPLING 
C RHO  - CORRELATION  COEFFICIENT  AT  EFFECTIVE  SAMPLING  RATE 

C SIGH Ah  - STANDARD  OEVIATICN  OF  NHITE-MOISE  I NFUT 

C INITIALIZE  GENERATOR 

IF  (I  .GT.  1)  CO  TO  1 
CALL  RANOOMRANNUH.l) 

XNOISt  s SaGHA*RANNUN 

C FIND  SAMPLING  TINE  TO  SATISFY  SHANNON  SAILING  THtORsN 
IF ((OELTAT/ TAU)  .GT.  1CI.)  GO  TO  J 
N = 1.*1C.»0£LTAT/TAU 
rho  * exp(>oelTat/(Tau*float {<«))) 

SIGHAh  * S:GHA*iQRTIl.-hHO**Z) 

RETURN 

3 CONTINUE 
« * l 
RHO  * 5. 

SIGHAh  s SIGH A 
RETURN 

1 CONTINUE 

Ou  Z J * ltN 

CALL  RA NOO" (RANNUNf I ) 

XNOISE  * RHC*  XNOlSf.  ♦SIGf  Ah*RANNUN 

2 CONI I HUE 
RETURN 
END 


C'ORKAL  - This  subroutine  is  the  general  (matrix)  version  of  the  Kalman  Filter  which  has 
been  rcstructed  to  accept  serially  correlated  measurement  error.  The  model  and 
algorithm  development  are  presented  in  Section  VI.  The  matrix  equations  which 
this  algorithm  computes  are  shown  in  Table  6.2. 


SUJEOUTINC  C0^<6Lt  *,PHlfFfP,Q,3,Hf  HnL0,PSI.Z,£T4,K,r«,.l,I0l;T) 

C LAST  UPDATE  - 1 NJVFHER  19X3 

C Gr NE°AL  <AlhAN  FTlT£R  CVClE  FOP  LInpAR  SYS  IF*  HIT*  AUTOCORRllATEC 
C MFASUPENENT  NO I St  - GIVEN  INPUT  AT  TTn£C«-1>,  THIS  POuTTwF  C£TU°NS 

C OUTPUT  »T  TIKOO  - MATRICES  A R£  N CT  RtDEFlNEO  EXCEPT  WHERE 

C INDICATED 

C XIN)  - STATE  VECTOR  - TNPUT«-1»  - 0UTPl|T(K| 

C PHI  IN, N)  - T = ANSriCN-  MATRIX  - INPUT (<,<- 1) 

C FIN)  - OCTt ■»«! NI  »T IC  cOsCING  VECTOR  - INFUT«-tl 
C P(*I»N>  - fRROP  COVAPTANCt  **ATRIX  - INPUTK-l/F-l)  - CUTSUT«/K» 

C OIN,N)  - PROCESS  NOISE  "ATPIX  - IN*UT(K-11 

C PIN, Hi  - *£ ASuRCNENT  pRRC*  MATRIX  - WH  IU  NOIsE  POKTTOM  - INPUT(K-I) 

C M(N,N)  - C-3j£AVATTON  MATRIX  - CUR^NT  - INPgTIO 
C M3LD(>*,N1  - OCSf’VATTON  MATRIX  - i.AST  T IH£  - NCX-11  - INPUTIX-U 

c - autocorrelation  ipansitton  matrix  - ikput<k,<-ii 

C Z(M)  - MC ASUR£HrN-  VcrT  OP  - INPUT {XI 

C ETAtHJ  - P=U1UCT  OP  =sicx.<-ll  tines  TK-  a POST£hIO®I  RESIDUAL  - 
C INPUT  K-l)  - OUT3UT(X) 

: - cits  -Am*  - calculated  - output i<) 

C « - SIZE  0r  X 
C *1  - SIZE  OP  Z 

c layr  - output  swit-h  iset  to  i to  dumr  hl  internal  output) 

C AlN»,C*(N,M),C(M,“>,TCN,M)  - W0P<INC  ?PACF 

C NjTt  - IN, HI  MUST  NUT  EXCEPT  (13, T*  - Jr  EJThLR  OOPS,  CHANGE  Tnr 
C DIMENSIONS  THE  WORKING  SP4CP  AN  C T Hu  *E  XT- *0 -t  A$T  ARGUMENT 

c op  **atk  jx  :*.vep$:cn 

•»£Al  x 

UIHp  N IION  » INI  , 0tiI  (N,N)  ,P(NI  ,P(»*,N>  ,Q(.N,N>  ,*  I M,  1*1  , H (H,N1  , 

Z ‘•vLUf.NI  ,PSI  If*,"!)  , Z(N>  , £T  A (**)  ,X(:J,M) 

OI-ENEIln  Alt,)  ,R(  if  ,1.|  ,C<  3,31  ,r(  P,3) 

C EXT p.APjl A 'ION 
C As  C*I »,»F 

C i s f «I  »o»  T9,iNScTSP  I phI  I »0 
tr  i I * :,s 
All)  * p ( ! I 
OP  1 J * 1,N 

All!  s A . 1 1 »OhT I !, J) *X I J» 

"II.Jl  - C d,  J) 

00  1 <C<  = 1 ,N 

00  1 L = 1 ,N 

911, J>  = u*I,J>  »PHI  (I,«1*P(«,L  1’PHllJ.L) 

1 CENTIME 

C C = P*h*=*  RANSPJSP'hi  »PSl»M0L0*c*TOANSA0S£IH0k.T)  * TRANSPOSE  1PSI1 

c -h»pki *p*t pan; pose  i “oloi*t®ansbo<e  :psd 

C -PS  I •f’C'LC*P*  TRANSPOSE  IPMl)  *T»S  l$P0SE  (HI 

u c > 1*1,- 

00  2 j - 1 ,- 

CU,J)  * P(I,I) 

00  Z <*  = 1,N 

oo  a l * i ,n 
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CCI,J>  = tK<)*8(ICK»l.  )*H(  «tL> 

00  ? N I * 1 
00  9 NN  s i,* 

0C  13  * i,<* 

IC^ITI  »U£*  * JLn*^ ■»*<*«>  *P  (*<,(.)•  PHI  (#*,  u*HOL  Of  J.HK) 

2 cr.HTiwr 

r 1 s_tJ*T,1N5‘*0Sf  »'O*®'«r»P«TRANSP0';E(M0Lr»*PSI 
00  M * Sit 
"C  3 J * 1,** 

tm,j>  s c. 

CO  3 « * 1,N 

0(1, Jl  s 0(I,J>  »!*(  It*0»H<  J <*) 

00  3 L s :,H 
00  3 in  S 1,« 

3 cc£rr «bc 

C 

: : * r<v»$f(o 

rr  I*  *r«*  ?*  3iL*  *C,4,0CT» 

, I^.f  ‘£0*  *•  CU,I)  * i./rti,o 

C < s y*C 

o>j  - I = j ,*i 

cc  » j * i,« 

«C(J,J)  s C. 

Ow  * <<  = i,h 

KlIIjJI  = X (I,J>*D(I,K<I*C(<«C,  J) 

<•  CCnTIMjf 
C '^OiTF 

0 < - S»<*  (2-m»  *«■£*■  » j 

00  5 I * i ,'j 
Ml)  * HH 
00  5 J s 1 ,•» 

*(tl  * *( I >-< (I , J) • (?( JJ ,£T* ( J) J 

CC  5 <«  * i,-j 

*«ij  = »ih*<ii,jj*ihj,<ki«iki(| 

3 CONTINUE 
C - = 5-3»< 

UC  b t s t,»l 
•JC  $ J * 1,N 
p<r,Ji  r 3 ( I , J) 

00  6 K<  a J,’* 

P»I,j)  s P(l,j).oir,r<»*i?cj,<<j 
6 CO*iTtPje 
C ;TS  r PSI*<£-H*<I 
OC  7 I * J,* 

rum  s t. 
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